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ABSTRACT 


This  thesis  is  a  collection  of  three  papers 
(Chapters  1-3)  and  a  research  note  (Chapter  4)  which  as 
of  November  16,  1979  have  the  following  dispositions. 
Chapter  one  has  been  published  in  its  entirety,  appearing 
in  the  June,  1979  issue  of  the  Bulletin  of  the  Seismolo- 
gical  Society  of  America;  chapter  two  is  forthcoming  in 
a  slightly  abridged  form  in  a  fall  issue  of  the  Canadian 
Journal  of  Earth  Sciences;  chapter  three  has  been 
submitted  for  publication  to  Geophysics;  and  chapter 
four  is  as  yet  unsubmitted.  Each  of  the  chapters  is  self- 
contained  and  does  not  require  the  presence  of  the  others. 

No  apology  is  made  for  the  treating  of  the  seemingly 
trivial  case  of  SH  waves  when  dealing  with  the  problems 
discussed,  as  it  is  the  author's  intent  to  present  in  the 
simplest  possible  form  certain  topics  in  elastic  wave 
propagation.  All  results  derived  here  can  and  in  some 
cases  have  been  applied  by  the  authors  to  more  sophis¬ 
ticated  topics  in  elastic  wave  propagation.  A  large 
part  of  this  thesis  may  appear  tutorial  in  nature,  but 
this  is  necessary  to  prepare  a  basis  from  which  to 
explore  obvious  extensions. 

There  are  many  reports  in  the  literature  of 
anomalies  in  travel  time  data  when  an  isotropic  homo¬ 
geneous  earth  model  is  used  to  interpret  field  data.  In 
several  instances,  the  introduction  of  a  layered  trans- 
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versely  isotropic  model  has  successfully  explained  these 
kinematic  irregularities.  However,  it  is  useful,  in  fact 
essential,  to  confirm  the  kinematic  fit  with  a  dynamic 
(amplitude)  comparison. 

An  asymptotic  expansion  of  the  displacement  vector 
in  terms  of  inverse  powers  of  frequency  is  employed  in 
chapter  one  to  investigate  the  dynamic  properties  of 
SH  waves  propagating  in  plane  layered  transversely 
isotropic  media.  Both  reflected  and  head  waves  are 
considered  in  terms  of  the  asymptotic  expansion,  and 
their  ranges  of  validity  are  discussed. 

Although  solutions  to  the  most  general  case  of  wave 
propagation  in  anisotropic  media  can  be  found  in  the 
literature  it  is  instructive  to  consider  this  simple  case 
in  which  many  quantities  inherent  to  wave  propagation  in 
anisotropic  media  can  be  more  readily  understood  and  can 
be  solved  analytically  rather  than  reverting  to  numerical 
methods. 

In  chapter  two  the  problem  of  SH  waves  propagating 
in  a  transversely  isotropic  plane  layered  medium  is 
discussed  through  the  use  of  integral  transforms  and 
evaluation  by  steepest  descents.  This  procedure  yields 
not  only  the  asymptotic  solution  which  is  also  attainable 
using  an  asymptotic  ray  series  approach,  but  also  allows 
for  the  investigation  of  the  interference  of  the  reflected 
and  head  waves  in  the  vicinity  of  the  critical  point 
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(point  of  critical  refraction) .  It  is  in  this  region 
that  asymptotic  ray  theory  breaks  down  or  in  the  least 
introduces  considerable  error  in  the  displacement 
amplitudes. 

The  theory  for  the  simple  case  of  SH  rays  propagat¬ 
ing  in  plane  layered  media  consisting  of  thick  layers 
separated  by  thin  layered  transition  zones  is  examined 
in  chapter  three.  The  SH  case  was  chosen,  as  the  basic 
idea  of  the  method  is  conveyed  without  an  excess  of 
mathematics  that  is  necessitated  by  the  introduction  of 
potentials  in  the  P-SV  case. 

A  stationary  phase  approximation  is  employed  when 
evaluating  the  integral  which  yields  the  displacement 
due  to  an  arbitrary  ray  propagating  in  the  above 
mentioned  medium,  and  the  validity  of  this  approximation 
is  discussed. 

The  comparison  of  ray,  numerical  integration  and 
ray-reflectivity  synthetic  sections  indicates  that  the 
transition  zone  method  gives  quite  acceptable  results  for 
small  source  receiver  offsets.  This  method  is  suitable 
for  application  in  the  oil  industry  as  individual 
arrivals  associated  with  ray  paths  in  the  thick  layers 
are  identified  in  the  synthetic  trace.  Furthermore, 
compared  to  other  methods,  this  technique  is  quite  cost 
efficient. 

In  several  papers  in  seismological  literature  it 


vi 


' 


' 

. 


. 


has  been  shown  mathematically  that  a  periodic  isotropic 
two  layered  medium  behaves  as  a  homogeneous  transversely 
isotropic  medium.  This  assumption  has  been  shown  to  be 
valid  as  far  as  the  propagation  of  elastic  waves  is 
concerned  provided  that  the  seismic  wavelengths  used 
are  large  when  compared  to  the  thicknesses  of  the 
individual  isotropic  layers.  Chapter  four  contains 
seismograms  computed  using  two  different  methods  which 
compare  the  kinematic  and  dynamic  properties  of  the 
periodic  structure  with  its  equivalent  isotropic  homo¬ 
geneous  structure  for  SH  waves.  Although  this  comparison 
seems  quite  trivial,  none  has  appeared  in  the  literature 
which  is  probably  due  to  the  prohibitive  computing  costs 
involved . 
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CHAPTER  1 


SH  WAVES  IN  LAYERED  TRANSVERSELY  ISOTROPIC  MEDIA 


AN  ASYMPTOTIC  EXPANSION  APPROACH 


1.1  Introduction 

An  asymptotic  expansion  of  the  displacement  vector 
in  terms  of  inverse  powers  of  frequency  is  employed  to 
investigate  the  dynamic  (amplitude)  properties  of  SH 
waves  propagating  in  plane  layered  transversely  isotropic 
media.  Both  reflected  and  head  waves  are  considered  in 
terms  of  the  asymptotic  expansion,  and  their  ranges  of 
validity  and  accuracy  are  discussed. 

Although  the  solution  of  the  most  general  case  of 
wave  propagation  in  an  anisotropic  medium  has  been 
presented  in  the  literature  (Cerveny  (1972)),  it  is 
instructive  to  consider  this  simple  case  in  which  many 
quantities  inherent  to  wave  propagation  in  anisotropic 
media  can  be  more  readily  understood  and  can  be  solved 
for  analytically  rather  than  reverting  to  numerical 
methods . 

1.2  Mathematical  Preliminaries 

In  a  transversely  isotropic  inhomogeneous  elastic 
medium,  neglecting  body  forces,  the  equation  of  motion 
for  the  propagation  of  an  SH  wave  can  be  written  as 


(1.1) 


' 


' 
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where  (x^x^x^)  are  the  rectangular  Cartesian  coordi¬ 
nates  and  ,  0^^  and  p  are  the  elastic  parameters  and 


density  required  to  specify  the  medium  for  SH  wave 
propagation  (Potsma  (1955),  and  Sato  and  Lapwood  (1968)). 
These  quantities  may  be  arbitrary  continuous  functions 
of  position. 

As  the  SH  displacement  vector  is  normal  to  the 
plane  of  incidence  in  a  transversely  isotropic  medium, 
the  displacement  vector  u  can  be  chosen  without  loss 
of  generality  as 

(1.2)  u  =  (0,u,0) . 

Although  the  problem  as  specified  by  equations  (1.1)  and 

(1.2)  can  be  treated  in  two  dimensions  (the  (x^,x^) 
plane) ,  it  is  mathematically  convenient  for  what  follows 
to  retain  the  derivative  with  respect  to  x2* 

Analogous  to  methods  used  in  geometrical  optics 
(Luneberg  (1944)  ,  Kline  (1951) )  the  solution  of  the 
equation  of  motion  will  be  assumed  to  have  the  form  of 
an  asymptotic  series  in  inverse  powers  of  frequency 
(Babich  and  Aleksiev  (1958)  ,  Karal  and  Keller  (1959)  , 
Cerveny  (1972),  Hron  and  Kanasewich  (1971))  so  that 


(1.3)  u(x1,x2,x3,t)  =  l 

n=0 


oo 


n  1  2  3 


n  i'  2'  3  SH 
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where 


oo 

-> 

n 


SH 


-  frequency 


-  a  unit  vector  in  the  direction 


t  (x^x^x^)  -  a  phase  function  which  describes 

the  wavefront  propagation 
Un  (xi ,xn 7X3)  -  complex  amplitude  which  is  only 

coordinate  dependent. 

The  substitution  of  the  assumed  asymptotic  series 
solution  (1.3)  into  the  equation  of  motion  (1.1),  yields 
with  the  added  conditions  that  U  ^  and  U  ^  identically 
equal  to  zero,  the  following  relation 


(1.4)  N  (U  )  -  M  (U  ,)  +  L  (U  9)  =  0 
n  n-1  n-2 


where 
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N(Ur)  (A44Pq  +  A44P2  +  A66P3  1)Un 


(1.5) 
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with  A 


44 


— —  ,  A  =  — —  and  p.  =  4^-  being  the 
p  66  p  1  ox. 


i-th  component  of  the  slowness  vector.  Discussions  of 
the  slowness  vector  and  slowness  surfaces  can  be  found 


in  Musgrave  (1970)  and  Kraut  (1963). 

Equation  (1.4)  must  be  equal  to  zero  for  all  values 
of  n  (n  =  O,00).  Thus  for  n  =  0,  remembering  that  U_1 
and  U_2  are  assumed  equal  to  zero,  the  following 
condition  must  hold 

(1.6)  ^44^1  +  ^44^2  +  "^44^3  ~  = 

Assuming  that  Uo  is  not  identically  zero  then 

(1.7)  G(xi,pi)  =  A44px2  +  A44p22  +  A66p32  =  1.  i  =  1,2, 

The  quantities  p^  are  the  slowness  vector  components 
and  are  related  to  the  wavefront  normal  velocity  (phase 
velocity)  through  the  relation  (Cerveny  (1972) ,  Cerveny 
and  Psencik  (1972) ) 

(.1.81  p.  =  ^  i  =  1,2,3 

where  Kb  are  the  directional  cosines  of  the  wavefront 
normal;  in  case, 

=  sinQcos^ 

(1.1)  N2  =  sin0sin4> 

=  cos0  . 

The  angle  0  is  the  angle  the  wavefront  normal  makes 
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with  the  x^  axis  and  <f>  is  the  angle  of  azimuth  of  the 
wavefront  normal  measured  from  the  x^  axis.  Thus 
utilizing  equations  (1.7),  (1.8)  and  (1.9)  the  wave- 

front  normal  velocity  can  be  written  as 

(1.10)  =  A-.sin^O  + A, rcos^0 

44  66 

from  which  it  is  obvious  that  the  wavefront  is  rota- 
tionally  invariant  about  the  x^  axis. 

1.3  Rays  and  Snell's  Law 

The  quantity  G(x^,p^)  defined  in  the  preceding 
section  is  a  homogeneous  function  of  order  2  in  p.. 
From  Eulers  theorem  on  homogeneous  functions,  the 
relations 

3 

(1.11)  l  p.  =  2G  can  be  obtained. 

i=l  1  3pi 

The  equation  G(x^,p^)  =  1  is  a  nonlinear  partial 
differential  equation  for  x(x^)  (which  describes  the 
propagation  of  the  wavefront)  and  may  be  solved  using 
the  method  of  characteristics  (Courant  and  Hilbert 
(1962)).  The  equations  of  the  characteristics  corres¬ 
ponding  to  the  partial  differential  equation 
G(x^,p^)  =  1  can  be  written  with  the  help  of  relation 
Cl. 11)  as 
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d.i2)  \ 


dx^ 

1 

3G 

A  .  .  sin6cos<j) 

-  A  r.  -  44 

dx 

2 

dp1 

> 

1 

i — 1 

h 

C 

1 

dx2 

1 

dG 

A  .  .  sin0coscj) 

—  A  T)  —  ^  ^ 

dx 

2 

dp2 

44P2  V 

dx3 

1 

dG 

A^  ^cos  0 
-An-  66 

dx 

2 

3p3 

66P3  V 

dx  . 

The  terms  are  interpreted  physically  as  the 

components  of  the  ray  velocity;  that  is,  the  velocity  of 
energy  propagation.  In  general,  in  an  anisotropic  medium, 
the  direction  of  energy  propagation  does  not  coincide 
with  the  wavefront  normal  direction  (Musgrave  (1970) ) . 

If  a  is  denoted  as  the  angle  the  ray  makes  with  the  x^ 
axis  (which  will  henceforth  be  assumed  to  be  the  vertical 
axis)  and  3  the  angle  the  ray  makes  with  the  x^  axis,  the 
relations  between  the  ray  angles  and  wavefront  normal 
angles  can  be  found  by  taking  the  ratios  of  the  compo¬ 
nents  of  the  ray  velocity  (equations  (1.12))  so  that 


tana  = 


L44 


A 


tan0 


66 


(1.13) 


tan3  =  tarup 

The  magnitude  of  the  ray  velocity,  a,  is  obtained 

dx 


from  the  inner  product  of  ^ 
itself ,  viz . , 


dx^  dx2  dx^ 


dx  '  dx  7  dx 


with 
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(1.14) 


a  = 


dx  _  dx 
dx  dx 


2  »  \  ^ 


A  ^sin^G  Ar  ^cos^0 
4  4 _  6  6 _ _  : 
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or  using  equation  (1.13) 


.  2  2 

, ,  ,  c  x  1  sin  a  .  cos  a 

(1-15)  ~2  -  ~a77~  +  ~^7T 

a  44  66 


Equation  (1.15)  is  that  of  an  ellipse,  signifying 
that  the  SH  wavefront  is  an  ellipsoid.  The  absence  of 
the  azimuthal  angle  6  from  (1.15)  implies  that  the  SH 
wavefront  is  an  ellipsoid  of  revolution,  the  x^  axis 
being  the  axis  of  symmetry.  It  should  be  noted  that  at 


G  =  a=0,  V  =  a=  /A66  and  at  6  =  a  =  |,  V  =  a  = 


44 


Two  new  variables  a  and  a  corresponding  respectively 

Z  X 


to  y^A66  and  yA^  ,  the  velocities  along  the  minor  and 
major  axis  of  the  ellipsoid,  will  be  introduced  so  that 


2  2  2  2  2 
(1.16)  V  =  a  sin  0  +  a  cos  0 

x  z 


(1.17) 


2  2 
sin  a  cos  a 


The  designation  of  the  vertical  axis  of  the  ellipsoid  as 
the  minor  and  the  horizontal  axis  as  the  major  is  done 
as  in  most  seismological  applications  a  >a  ;  that  is. 


x  z 
a 

x 


the  degree  of  anisotropy  of  the  medium  — 
greater  than  one. 


is  usually 
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The  geometry  of  a  part  of  an  SH  wavefront  in  an 
inhomogeneous  transversely  isotropic  medium  is  shown  in 
Figure  1.  It  can  be  seen  that  the  projection  of  the  ray 

velocity  onto  the  wavefront  normal  yields  the  magnitude 

of  the  normal  velocity,  so  that 

(1.18)  V-a  =  v2. 


This  same  result  can  be  obtained  using 


V 


and 


-v 

a 


a 

x 


) 


together  with  equations  (1.8)  and  (1.9).  By  definition 
of  the  inner  product,  V«a  =  Va  cos (0-a)  so  that 

(1.19)  a  cos(0-a)  =  V. 

1.4  Calculation  of  Ray  Amplitudes 

The  solution  for  the  complex  amplitude  components 

Un^xl'x2'x3^  ^he  raY  series  expansion  (1.3)  is  an 

iterative  one;  that  is,  if  the  solution  for  U  is  to  be 

found  all  U^,  k<n  must  be  known.  In  equation  (1.4)  it 

can  then  be  assumed  that  L(U  , )  is  known.  Call  this 

n-1 

2  2 

quantity  qn*  From  (1.4)  with  +  ^44^2  + 

2 

A^p-,  -  1)  =  0  the  following  expression  results: 

b  b  J 


' 

Geometry  of  the  wave  front  and  ray  and  the 
relation  between  the  wave  front  normal 
velocity  and  the  ray  velocity  a. 


Figure  1.1 
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3U 


2  A  .p 


n 


3U 


4  4  1  3x 


+  2A .  p 


n 


1 


442  3x. 


+  2A  p 


3U  U 
n  +  n 


663  3x3  p 


(1.20) 


3^  (pA44Pl)+dq  (pA44P2)  +  3^ 


=  q 


n 


With  the  aid  of  equations  (1.12)  this  can  be  written 


I  3U  dxn  3U  dx_  3U  dx_ 
2  n  1  +  n  2  +  n  3 


[_3x,  3t  3x„  di 


(1.21) 


3x^  dx 


U 


n 


i  3x 


3  (pai>  +  ?lr  (pa2>  +  db  <pas) 


n 


where  a^f i  =  1,2,3  are  the  components  of  the  ray  velocity. 
Equation  (1.21)  can  be  written  more  compactly  as 


dU  U  q 

(1-22)  dF  +  2?  7-<pa»  -  f 


Employing  a  vector  integral  theorem  due  to  Gauss 
as  demonstrated  in  Appendix  A1  equation  (1.22)  becomes 
the  transport  equation 


dU  u  q 

(1-23)  dF  +  im  If  <pVJ>  -  -f 


the  quantity  J  being  defined  in  the  appendix.  The 
solution  of  (1.22)  is  given  by 


(1.24)  U  ( t) =U  (t  ) 
n  n  o 


(pvj)  t 
Lo 


(PVU) 


+ 


t  J  2  (  PVJ)  t 


(pVJ)T2  qn(T)dx 


■ 


' 


11 


where  t  defines  a  reference  wave  surface,  at  which  the 
o  ' 

value  of  U  (t  )  is  an  initial  condition  which  is  assumed 
n  o 

known.  The  integration  is  taken  to  be  along  the  ray  path 
from  the  reference  wave  surface  to  an  arbitrary  wave 
surface  denoted  by  t. 

Consequently  from  formula  (1.24),  Un(t)  can  be 
determined  at  an  arbitrary  point  along  the  ray  if  its 
value  at  a  previous  point  along  the  ray  is  known.  In 
seismological  applications  the  zero  order  or  leading 
term  in  the  ray  expansion  is  usually  considered  to  be 
a  sufficiently  good  approximation.  As  qQ  =  0,  equation 
(1.24)  becomes  for  n  =  0 


(1.25)  UQ(t)  =  UQ(to) 


(PVJ) 


(pVJ) 


In  the  problem  under  consideration  p,  V  and  J  are  all 
easily  computed  along  the  ray. 

Inspection  of  equation  (1.25)  reveals  that  an 
infinite  value  of  the  ray  amplitude  is  predicted  at 
caustics,  where  the  function  J,  which  is  related  to  the 
Jacobian  of  the  transformation  from  Cartesian  coordinates 
(xlfx2,x3)  to  ray  coordinates  (a,3,x),  vanishes. 

Although  some  increase  in  amplitude  in  the  vicinity 
of  a  caustic  may  be  physically  justified  by  claiming 
that  energy  flux  becomes  highly  concentrated  when  a  ray 
tube  collapses  to  a  condition  of  zero  cross  sectional 


' 


■ 


12 

area,  the  inapplicability  of  equation  (1.25)  in  the  region 
of  a  caustic  steins  from  purely  mathematical  considerations. 

It  is  well  known,  for  example,  Babich  and  Buldirev 
(1972)  or  Cerveny  et.  al.  (1977),  that  the  ray  series 
expansion  in  equation  (1.3)  can  be  used  only  in  such 
regions  where  the  ray  field  is  regular,  implying  a  non¬ 
zero  value  of  the  Jacobian  D (x^ ,x2 ,x^) /D (a , 3 , t ) .  As  this 
necessary  condition  for  the  application  of  the  ray  series 
expansion  technique  is  clearly  violated  at  a  caustic,  no 
attempt  should  be  made  to  discredit  asymptotic  ray  theory 
by  extending  its  results  into  this  region. 

Fortunately,  alternative  high  frequency  techniques 
have  been  developed  for  the  evaluation  of  ray  amplitudes 
in  the  vicinity  of  a  caustic.  Numerous  papers  on  these 
techniques  appear  in  the  literature  and  are  well  docu¬ 
mented,  for  example,  in  the  monographs  of  Stravroudis 
(1972)  ,  Babich  and  Buldirev  (1972) ,  Babich  and 
Kirpichnikova  (1974)  and  Cerveny  et.  al.  (1977). 

As  the  results  produced  by  these  alternative  methods 
link  smoothly  with  the  asymptotic  ray  theory  solution  in 
the  regions  where  the  latter  is  justified  (Choi  (1978) ) , 
equation  (1.25)  can  be  used  at  any  point  of  the  regular 
ray  field,  provided  that  the  effect  of  each  caustic 
which  the  ray  passes  through,  is  included.  In  the  high 
frequency  limit  this  effect  amounts  to  a  change  in 
phase  of  ±  2  for  each  caustic  passed  through. 


* 


' 

. 
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1.5  Reflection  and  Transmission  Coefficients 

At  an  interface  between  two  transversely  isotropic 
media  an  SH  wavefront  impinging  on  the  interface  from 
either  medium  gives  rise  to  a  reflected  and  transmitted 
wavefront.  Let  0^  (upper  medium) and  0~  (lower  medium) 
be  the  acute  angles  the  incident  SH  wavefront  normals 
make  with  the  vertical  axis  and  0^,  v  =  1,2  be  the  acute 
angles  the  reflected  and  transmitted  wavefronts  make  with 
the  vertical  axis.  The  reflected  or  transmitted  wave- 
front  in  the  upper  medium  will  be  specified  by  v  =  1 
while  v  =  2  will  typify  the  same  quantities  in  the  lower 
medium,  as  shown  in  Figure  2.  To  simplify  notation  in 
what  follows,  the  horizontal  axis  will  be  denoted  as  the 
x-axis  and  the  vertical  asix  as  the  z-axis  with  the 
positive  direction  chosen  to  be  downwards.  Further,  it 
will  be  assumed  that  the  interface  is  planar  and  that 
the  axes  of  rotational  invariance  of  the  wavefronts  are 
aligned  perpendicular  to  the  interface.  Thus  Snell's 
Law  for  this  case  can  be  written 


sin0 


Cl. 26) 


0 


sin0 


v 


V. 


V 


v 


where  it  is  to  be  remembered  that  Vq  and  are 
functions  of  0^  and  0^.  Employing  the  equations  (1.13), 
(1.14)  and  (1.15)  a  modified  Snell's  Law  can  be  obtained 
in  terms  of  the  ray  velocities  and  angles  in  the  form 


- 
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Figure  1.2  The  wave  front  normals  of  incidence  (0^  and 

0q)  and  reflection  and  refraction  (0^  and  0^) 
of  SH  wave  fronts  at  an  interface  between  two 
transversely  isotopic  media. 


15 


(1.27) 


aosinao 


a  sina 

v  v 


x 


0 


V 


If  the  two  media  described  above  are  in  welded 
contact,  the  two  boundary  conditions  which  must  be 
satisfied  at  their  common  interface  are  the  continuity 
of  shear  stress  and  displacement.  The  plane  of  inci¬ 
dence  can  be  chosen  without  any  loss  of  generality  as 
the  (x,z)  plane.  Under  the  last  assumption,  an  arbitrary 
displacement  vector  in  either  medium  can  be  written  as 

00  U  (x  ,  z )  exp  [  ico  ( t-T  )  ] 

(1.28)  u  ( x ,  z  ,  t )  =  l  - - - 

v  n=0  (iw) 


where  U  (x,z)  =  U  (x,z  )nptJ,  nOTT  being  a  unit  vector 
nv  nv  SH'  SH  3 

perpendicular  to  the  plane  of  incidence  (in  the  y 
direction) .  The  subscript  v  identifies  the  type  and 
medium  of  the  displacement  vector  while  does  the  same 
for  the  wavefront. 


v  =  0  -  incident  from  medium  1 

v  =  0  -  incident  from  medium  2 

v  =  1  -  reflected/transmitted  in  medium  1 

v  =  2  -  reflected/transmitted  in  medium  2 


Continuity  of  displacement  has 


r 


(1.29) 


, 
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It  should  be  recalled  that 


9  T 


(1.30) 


v 


sin0 


v 


9x 


V 


v 


and  that 


(1.31) 


9t  V+6  COS0 

v  _  vO  v 


9x 


V 


V 


The  choice  of  sign  in  (1.31)  is  determined  by  the  physical 
requirement  that  the  solution  for  the  disturbance  approach 
zero  as  the  wavefront  moves  away  from  the  interface  to 
infinity. 

The  continuity  of  shear  stress  requires 


(1.32) 


r 


l  (-1)V  o  (u  )  = 


a  (uj 
yz  0 


v=l 


yz  v ' 


.-°yz<U6) 


where 


a  (u  )  =  C  .  . 
yz  v  44 


(v) 


9 (u  )  a (u  ) 
y  v  z  v 


9  z 


9Y  _ 


This  yields 


(1.33)  l 
v=l 


2  C,  .  (V) cos6  U  , 
44  v  nv 


V 


v 


2  ,  s  9U  -  n  > 

F  +  I  (-1)V  C  (V)  (n'1)V 

v=l 


’44 


9  z 


where  for  incidence  from  above 
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C 44  cos6 

p  =  _lz - H  n  -  c 

VQ  u0  44 


(°)  aU(n-l) 0 
9  z 


and  for  incidence  from  below 


(0) 


F  = 


'44 


COS0 


V 


0 


0  (6)  3U(n-l)0 

nO  44  3z 


In  the  zero  order  approximation  (n=0)  equations 

(1.29)  and  (1.33)  yield  two  sets  of  linear  equations  in 

two  unknowns;  the  unknowns  being  the  ratios  of  the 

reflected  and  transmitted  amplitudes  to  the  incident 

amplitude.  At  this  point  it  should  be  recalled  that  in 

an  earlier  section  U,  ,  was  set  equal  to  zero. 

(-l)v  ^ 

As  the  axis  of  rotational  invariance  was  chosen 


perpendicular  to  the  interface,  it  follows  that  8q=0^, 
and  9q=02  an<^  consequently  vq=vj_  an<3  vg=V2* 


At  the  interface  C 

r  <°>  =  r  (2)  =  r 
C44  C44  B2* 


(0)  _ 


(1) 


44 


Thus 


(1.34a) 


-1 


B^cos6^  B2cos02 


V. 


V, 


'44 


U 


01 


U 


02 


=  B^  and 


-1 


B^cos0  ^ 

vT 


U 


00 


for  a  wavefront  incident  from  above  (medium  1)  and 


* 


' 
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(1.34b) 


1 

-1 

1 

1 — 1 

0 

D 

L 

1 

B^cosB^ 

B2COS02 

1 

a 

0 

NJ 

1 _ 

B2cos02 

1 — 1 
> 

_ l 

V2  j 

V2 

U 


00 


for  incidence  from  below  (medium  2) .  The  solutions  of 
(1.34a)  and  (1.34b)  are  given  as  R  where  y  typifies  the 

y  i 

incident  and  v  the  resultant  wavefront  moving  away  from 
the  interface. 


R 


B^cosB^  B^cosB^' 


11 


V. 


V. 


/D 


'12 


2B^cos0^ 

Vt 


/D 


(1.35)  ( 


'21 


2B2cos02 

VZ 


/D 


R 


22 


B  cos0.  B..COS0.,  \ 
Z  Z  1  1  \ 


V. 


V. 


/D 


D  = 


B^cosB^  B2cos02 


V. 


V, 


1.6  Geometrical  Spreading 

The  quantity  J  mentioned  in  equations  (1.23),  (1.24) 

and  (1.25)  can  be  physically  interpreted  as  the  measure 
of  incremental  cross  sectional  area  of  a  ray  tube  (as 
shown  in  the  appendix)  and  consequently  it  is  conceptually 
convenient  to  write  equation  (1.25)  in  the  form 
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P 


(1.36)  U (M)  =  U(M  ) 


V (Mq) p (Mq ) 1 \ 


V  (M)  p (M) 


1% 

da  (M  ) 
o 

1 

J 

_da  (M)  _ 

where  M  is  a  point  on  the  reference  surface  and  M  is  an 
o 

arbitrary  point  along  the  ray. 

If  the  ray  traverses  through  a  medium  from  a 

reference  point  M  and  encounters  an  interface  at  the 

o 

point  0  and  is  either  reflected  back  into  the  medium  of 
incidence  or  transmitted  into  the  second  medium,  the 
amplitude  measured  by  a  receiver  at  point  M  is  generally 
different  than  that  given  by  equation  (1.36)  as  the 
discontinuity  of  the  elastic  parameters  at  the  interface 
introduces  effects  which  must  be  taken  into  account.  The 
ray  now  essentially  consists  of  two  segments  and  the 
amplitude  at  M  can  be  written  as 


(1.37) 


U  (M)  =  U(M  ) 
o 


V(M  )  p  (M  )~h  dp  (M  )1% 
o  o  !  o 

V  (M)  p  (M )  j  _dp(M)  J 


R 


where  V  and  p  denote  the  velocity  and  density  at  the 
point  0  on  the  side  of  the  incident  ray,  V1  and  p'  denote 
the  same  parameters  at  the  beginning  of  the  second  seg¬ 
ment  of  the  ray  and  R  is  the  appropriate  reflection  or 
transmission  coefficient. 

Assuming  that  a  total  ray  is  composed  of  k+1 


* 
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segments,  where  a  ray  segment  will  be  defined  as  the  part 
of  the  total  ray  between  one  reflection  or  refraction  and 
the  next  (Figure  3),  equation  (1.37)  can  be  generalized 
to  read 


(1. 38) 


U(M) 


U  (M  ) 
o 

L 


V  (M  )  p  (M  r  ^  k 

— 9. - 2_!  n 

V  (M)  p  (M)  J 


V'  p 
Vp 


R  . 
3 


where  the  primed  quantities  are  as  previously  defined, 
Rj  is  the  reflection  or  transmission  coefficient  at  the 
point  0^  and  L  is  given  by 


(1.39) 


L  = 


da(M)  | 

da  (M  )  j 

oj 


>  k 

n 

j=i 


da  2 

L^Jo. 


and  serves  as  a  measure  of  the  geometrical  spreading  of 
the  ray  tube. 

The  remainder  of  this  section  will  be  devoted  to 
obtaining  an  expression  for  the  geometrical  spreading  L, 
for  the  simple  case  of  plane  parallel  homogeneous  layers, 
in  which  the  axes  of  rotational  invariance  of  the  wave- 
fronts  are  perpendicular  to  the  interfaces. 

If  a  and  3  are  parameters  which  describe  the  ray, 
da  can  be  expressed  using  the  standard  formula  from  the 
differential  geometry  of  surfaces  as 

(1.40)  da  =  J  dad3 


where 


' 

g  f 


' 


Figure  1.3  Definition  of  interface  coordinates  and  ray 


segments . 
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2  h 

J  =  (EH  -  G  )  ‘ 

8R.3R  r  =  8^.8$  8R  9R 

h  da'doc  b  3ct*3B  H  3  3*93 

”)■  ,  t 

R  =  the  position  vector  of  the  ray  of 
magnitude  R. 

As  the  SH  wavefront  is  an  ellipsoid  of  revolution  in 

0 

a  homogeneous  medium,  the  components  (x,y,z)  of  the  vector 
R  can  be  expressed  in  standard  polar  coordinate  notation 
as 


(  x  =  RsinacosB  =  at  sinacosB  = 


a  t  sinacosB 
x _ 

{F+ ( 1-F) sin^a} 2 


(1.41)  y  = 


v  z  = 


a  t  sinasinB 
x _ 

{F+  (1-F) sin^a}^ 


a  t  cosa 
_ x _ 

{F+ (1-F) sin2a}^ 


2 

where  a  =  a  /{F+ (1-F) sin  a} 2  is  the  ray  velocity,  a  is 

X  X 

the  velocity  of  the  wavefront  in  the  horizontal  direction 

((x,y)  plane),  a  is  the  velocity  of  the  wavefront  in 

z 

2 

the  vertical  (z)  direction  and  F  =  (a  /a  )  . 

X  z 

Implementing  equations  (1.40) 

2  2  2  2  2 
a  t  (sin  a+F  cos  a) 

h  “  2  3 

(Ft(l-F)sin  a) 


G 


0 


' 
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H 


2  2  .  2 
a  t  sm  a 
x 


(F+ ( 1-F) sin^a) 


and  thus 


(1.42) 


da 


2 

R  smadadg 
cos(e-a)  ' 


0  being  the  angle  which  the  wavefront  normal  makes  with 
the  vertical  axis. 

The  reference  surface  will  be  chosen  such  that  the 
component  of  the  position  vector  defining  a  point  on  the 
wave  surface  has  unit  length  in  the  vertical  direction. 
This  being  so,  with  R  =  at,  t  must  have  a  magnitude  of 

— - — .  Thus,  R  =  — - — ,  where  a  has  units  of  t  \ 
i  a  |  o  |  a  |  z 

|  z  |  I  z  I 

This  yields 

2 

a  sina  da  df3 

(1.43)  da  (M  )  =  - = — 2 — - — —  . 

O  I  I  Z  /  r\  \ 

a  cos  ( 0  -a  ) 

1  z  o  o 


An  alternate  expression  for  da(M)  can  be  obtained 
from  simple  geometrical  considerations  (Figure  4) . 

(1.44)  da (M)  =  d£1d£2 


where  it  can  be  seen  from  Figure  4 


(1.45) 


r 

dJL 

< 


cosa„da 
M  o 


cos  ( 0.  -a,.) 
M  M 


r  d3 

o 


■ 


Figure  1.4 


The  geometry  of  a  ray  tube  and  definition  of 
quantities  used  to  calculate  an  alternate 
expression  for  da(M)  [Equation  (1.45)]. 
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which  yields 


(1.46)  L 


1 


'1  r 


9r 


9a 


cosaM  cos (6  -a  ) 

_ M  O  O  ! 

sina  cos  (6.  -a..)  ! 
o  M  M  I 


k 

n 

j=i 


i  da 

LFj  0 . 

3 


where  it  has  been  assumed  that  M  lies  in  layer  1 
accounting  for  the  subscripts  on  a and  a  and  further, 
that  da (M)  is  valid  even  if  the  point  M  does  not  lie  in 
the  same  layer  as  Mq. 

It  can  be  seen  from  Figure  5  that 


Cl. 47) 


r 


-- 1 x. 


da 


_da  jo. 


cosa 

cosa' 

J  j 


cos  ( 6 1 -a  * )  ^ 
cos (0-a)  Q 


and  after  introducing  this  interface  correction, 
equation  (1.46)  becomes 

2 


(1.48)  L  = 


J1  r 


9r 


9a 


cosa^  cos(0o~ao) 
sinaQ  oos(0M-o^) 


k  r 
n 

j=i 


cosa 


L. 


2  cos  ( 0  1  -a  * ) 
cosa’  _  cos  (0-a) 

0i~  __0. 

9  9 


In  the  case  under  consideration,  the  interfaces  are 
plane  and  the  media  homogeneous  so  that  a'(CK  ^)  =  a(O^) 
and 


k 

n 

j=i 


h 

,  ,~]h 

cos ( 0  -a ' ) ; 

cosa 

1  cosa 

o 

j_cosa 1  Q 

cos (0-a) 

_°°SaM. 

(eM-v 


To  aid  in  numerical  computations  it  is  convenient 


to  parameterize  all  angles  in  terms  of  one  angle,  say 


■ 


- 


Figure  1.5  The  discontinuity  of  a  ray  tube  at  an  inter¬ 
face  indicating  the  interface  correction 
which  must  be  introduced  into  the  expression 
for  geometrical  spreading.  It  can  be  seen 
that  at  the  interface  da/da '  =  d£/d£ 1 ,  as  the 
spreading  in  the  azimuthal  direction  is  the 


same  on  both  sides  of  the  interface. 
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aQ  =  a^.  If  the  acute  angle  the  j-th  ray  segment 


makes 


with  the  vertical  is  denoted  as  ou  then  from  equation  (1.27) 


a  sina  a  .  sina  . 
(1.27)'  -± - ^-±-=-3 - =-1 


X- 


X  . 

3 


where  the  ray  velocity  along  the  j-th  ray  segment  is  given 
by 


(1.50) 


a  .  = 
3 


x  . 
1 


.a  \  2 

( 


(F  .  +  ( 1-F . ) sin  a . ] 
3  3  3 


— t - ,  F.  = 

.  2  j  i  a  I 


From  the  above  two  equations  the  following  quantities 
can  be  obtained 


sma.  = 


/  k . F . sina, 
^  33...  1 


3  [1+ (F.-l)k. sin2an 


(1.51)  s  cosa  . 

!  3 


n  .  •  2 
1-k.sin  a. 

3  1 


l 


k. 

3 


1+ (F . -1) k . sin  a, 
3  3  1 


2  2 
a,  a 
1  x . 

2 


For  plane  homogeneous  layers 


k+1 


k+1 


(1.52) 


=  I  rj  =  I 


j  =  l 


j  =  l 


h .  tana . 
3  3 


k+1  h.  /  k  .  F  .  sina.. 
y  3^  3.3..  1 

j  =  l  [  1-k . sin  a.  ) 

J  V  3  l; 


- 

where  hj  is  the  thickness  of  the  layer  in  which  the  j-th 

ray  segment  is  propagating  and  r_.  is  the  horizontal 

distance  the  j-th  ray  segment  travels. 

As  an  analytic  expression  for  r.  in  terms  of 

exists  the  calculation  of  -5 -  =  77 -  presents  no  problems 

9a  9a, 
o  1 

so  that 


(1.53) 


9r 

9a, 


2  2 

,,,  h.a  a/ 

k+1  3  x.  1 

j  =  l  a  a *  2a  2 

Zj  zx  xx 


=  I 


cosa. 


.  .  .  2  \3/2  * 

1-k  .sin  a, 

3  1 


Substituting  the  expressions  (1.49),  (1.52)  and  (1.53) 

into  equation  (1.48),  the  resultant  expression  measuring 
the  geometrical  spreading,  equation  (1.54)  becomes 


(1.54)  L  = 


~2  cosa^ 


x. 


h.a 
3  x 


1 


k+1 

T  I 


h.a 
3  x. 


Jl 


k+1 

I  o  57  L  o  0/9 

j=l  a  (1-k. sin  a,)  2  j=l  a  (1-k. sin  a,)  ' 

J  z .  3  1  J  z.  3  1 

3  J  3 


For  homogeneous  media,  the  density  and  normal 
velocity  is  constant  along  each  ray  segment.  This 
results  in  equation  (1.38),  the  expression  for  the 
complex  amplitude,  reducing  to 

U (M  )  k+1 

(1.55)  U  (M)  =  — rrS-  n  R. 

L  j=l  3 


where  L  is  given  by  (1.54). 
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1.7  The  Head  Wave 

In  this  section  the  head  wave  corresponding  to 
critically  refracted  ray  will  be  examined.  This  problem 
for  the  isotropic  case,  including  SH  waves  as  well  as 
coupled  P-SV  waves  has  been  treated  extensively  in  the 
literature,  notably  the  book  by  Cerveny  and  Ravindra 


(1971)  in  which  numerous  references  are  cited.  The 
approach  used  will  employ  asymptotic  ray  theory  which  as 
the  name  suggests  asymptotically  approaches  the  solution 
attainable  by  exact  (integral  transform)  methods.  This 
asymptotic  solution  is  reasonably  accurate  at  distances 
removed  from  the  critical  point,  however  it  tends  to 
break  down  in  the  vicinity  of  this  point.  Other  methods 
of  solution  must  be  used  here,  due  not  only  to  the 
inapplicability  of  asymptotic  ray  theory  but  also  the 
effects  introduced  due  to  the  interference  of  the 
reflected  and  head  wave  near  the  critical  point  (Cerveny, 
(1965) ,  Cerveny  (1967) ) . 

Considered  will  be  two  half  spaces  in  welded 
contact,  the  parameters  of  the  upper  and  lower  media 
being  denoted  1  and  2,  respectively.  The  velocity 
structure  will  be  assumed  such  that  for  some  real  angle 
of  incidence  a*,a^  =  ^  so  that  from  equation  (1.27) 


sina* 


1  a  2 
X1 


which  is  equivalent  to  saying  that  for 


7T 

0*  (wavefront  normal  angle)  ®2~2 


some 


so  that 


. 

* 
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sin  6^* 


A  reference  point  M  will  be  chosen  a  distance  h 

o 

above  the  interface  at  t=0  and  a  two  segment  ray,  the 
geometry  of  which  is  shown  in  Figure  6  will  be 
considered.  Then 


and 


r  =  h  tana,  +  H  tana,, 
_L  A 


(1.56) 


£  =  -  =  (r  -  h  tana,  )/sina, 

cosa„  1"  2 


3  IT 

(1.57) 


2 

cos  a. 


2  2  3 

a  a  a,cosan 

x2  Z2  1  1 

2  2  3 

a  a  a 
x^  z^  2 


£  + 


,  2  2  3  2 

ha  a  a,cos  a, 
x1  z1  2  2 

2  2  3  3 

a  a  a, cos  a, 
x2  z2  1  1 
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2 

cos  a2 


3r 
3  a^ 


so  that  the  measure  of  the  geometrical  spreading  from 
equation  (1.48)  is 


(1.58) 


L  = 


1 


cosa, 


1  r 


3r 

3a, 


etna. 


cosa. 


The  expression  for  the  zero  order  term  approximation 


Figure  1.6  The  geometry  and  notation  used  in  deriving  an 
expression  for  the  head  wave  amplitude.  The 
points  A+  and  A-  lie,  respectively,  above  and 


below  the  interface.  To  obtain 
at  A-  in  terms  of  the  amplitude 
transmission  effects  across  the 
must  be  taken  into  account. 


the  amplitude 
at  A+ ,  the 
interface 


. 
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of  the  complex  amplitude  at  the  point  A  is 

R12 

(1.59)  UQ(A)  =  — “  cosa2  =  Uq(A)  cosc<2 

L 

where  R^2  is  the  transmission  coefficient  given  by 
equation  (1.35).  From  (1.59)  it  is  obvious  that  for 

7T 

a2  =  2  (critical  refraction) ,  the  complex  amplitude  at 
the  point  A+  just  below  the  interface  is  zero,  and  as  a 
consequence  to  determine  the  amplitude  at  A+ ,  higher 
order  terms  in  the  ray  series  expansion  must  be 
considered . 

It  will  be  shown  in  what  follows,  that  the  first 
non  vanishing  term  or  leading  term  in  the  series  is 

(.1-60)  u(A+)  =  YTwTU1^A+^  exp  ia)^  t t_T12  (A+)l  J  ^SH’ 

The  critically  refracted  ray  can  be  considered  as  a 
special  case  of  the  ray  incident  at  the  interface  from 
medium  2.  Before  proceeding,  it  is  useful  to  take  into 
account  the  interface  effects  and  obtain  an  expression 
for  this  critically  refracted  ray  at  A_  just  above  the 
interface  in  the  upper  half-space. 

To  determine  the  amplitude  just  above  the  interface 
at  A  in  medium  1  the  effects  of  transmission  across 
the  interface  from  medium  2  to  medium  1  must  be  deter¬ 
mined.  For  02  =  a2  =  \  the  cont:*-nuity  of  displacement 


►  ■ 


*  \  $■  ■'  ’ll 
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and  shear  stress  at  the  interface  for  incidence  from 
medium  2  using  equations  (1.29)  and  (1.33)  yields 


r 


-l 


(1.61)  ! 


B^cos6^  B^cosB^ 


L  v. 


V. 


w 

0 

1 

1 

1 

1 

< 

1 — 1 
D 

_ 1 

wv 

-  3  z  3 

ft*-77 

62"2 


where  it  is  to  be  remembered  that  6.  and  V.  (j=l,2) 

3  3 

refer  to  the  wavefront  and  not  the  ray. 

The  solution  of  (1.61)  for  the  complex  amplitude 
of  the  first  term  in  the  ray  series  in  the  ray  refracted 
back  into  medium  1  is 


B  au  (A  >  M  3u  (A  ) 

(1.62)  U,  (A  )  =  - - 5 - —  =  - - - — 

1  -  B^coseJ/V*  3z 


3  z 


Thus  the  expression  for  the  displacement  vector  at  a 
point  A  just  above  the  interface  is 


(1.63) 


M  3U  (A  ) 

u*  (A  )  =  ;  -x — ttz-  exp {  ioo  [  t- T 1 2 i  ( A_ )  ]  } n 


(iw)  3z 


SH' 


But  from  equations  (1.57),  (1.58)  and  (1.59) 


(1.64) 


3U0  <A+)  lim  3Uq  (A+)  x  lim  3U0  (A+) 
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This  ray  which  is  critically  refracted  at  the  inter¬ 
face  between  the  two  halfspaces  and  transmitted  back  into 
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the  first  halfspace  is  called  a  head  wave  (Cerveny  and 

Ravindra  (1971)).  Introducing  the  head  wave  coefficient 
* 

=  Ri2M21  t*ie  ^eac^  wave  amplitude  at  A  can  be 
written  as 


(1.65)  u* (A  )  = 


r  f 
121  2 


^tana* 


tor  (A  )  £ 


3/2  exp{ [ t-i121 (A_) ]  2>  ngH 


To  compute  the  displacement  of  the  head  wave  due  to 
the  leading  term  in  the  ray  series  at  the  point  M  in 
medium  1,  the  following  relation  can  be  used 


|-da(A  )  ^ 

(1.66)  UMM)  =  [_da(M)'_ 


As  the  waves  emanating  from  A_  are  conical,  geometrical 
spreading  occurs  only  in  the  azimuthal  (3)  direction  so 
that 

da(A_)  r(A_) 

(1*67)  da (M)  =  r (M) 


and  the  head  wave  displacement  vector  at  M  is  given  by 


(1.68)  u* (M)  = 


ri21F2^tana* 


2/2~  exp{ico  [t-i121  (M)  ]  - 


ITT  ,  -* 

2  }  nSH' 


wr  (M)  £ 
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1.8  Numerical  Results 

Figures  7  and  8  demonstrate  the  effects  that  varying 

anisotropy  has  on  reflection  and  transmission  coefficients. 

The  parameters  a  and  a  are  varied  in  the  upper  medium 

X  z 

(medium  1)  while  the  velocities  in  the  second  medium  are 

held  constant.  The  coefficients  are  plotted  for  0%,  10% 

and  20%  anisotropy  in  medium  1  and  10%  anisotropy  in 

medium  2.  The  degree  of  anisotropy  in  both  media  is 

/ a  -a  \ 

defined  to  be  ^ ^  x  100%.  The  actual  velocities 

a 

z 

used  for  each  of  the  three  curves  are  given  in  Table  1. 

An  attempt  has  been  made  to  negate  the  effects  of 

changing  impedance  at  the  interface  by  forcing  (a  +a  ) 

X1  Z1 

to  be  constant  in  each  of  the  three  cases  considered. 

As  the  coefficients  are  complex  quantities,  two  plots, 
one  of  ray  angle  of  incidence  vs.  amplitude  and  one  of 
ray  angle  of  incidence  vs.  phase,  are  presented  for 
each  of  the  four  coefficients.  The  phase  is  plotted 
only  once  (for  10%  anisotropy)  as  it  is  obvious  how  it 
will  behave  in  the  other  two  cases. 

A  comparison  of  log.  amplitude-distance  curves 
for  a  once  reflected  ray  and  its  corresponding  head 
wave  is  given  in  Figure  9.  The  curves  lettered  (a) 
correspond  to  0%  anisotropy  in  layer  1  while  those 
lettered  (b)  refer  to  20%  anisotropy  in  layer  1  and  the 
velocities  used  are  those  given  in  Table  1.  The  loga- 
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Figure  1 . 7 


The  reflection  coefficient  H1H1  and  the 


transmission  coefficient  H2H1  for  0(i), 
and  20(a)  percent  anisotropy  in  layer  1, 
for  the  media  specified  in  Table  1. 
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The  transmission  coefficient  H1H2  and 
reflection  coefficient  H2H2  for  0(i), 
and  20(a)  percent  anisotropy  in  layer 
the  media  specified  in  Table  1. 
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Figure  1 . 8 
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rithmic  singularities  in  the  reflected  amplitudes  are  due 
to  a  zero  in  the  reflection  coefficient  H1H1  as  shown  in 
Figure  7a. 

At  the  critical  points  (denoted  by  v's  in  Figure  9) 
the  head  wave  amplitude  blows  up  due  to  the  fact  that 
the  asymptotic  expansion  is  proportional  to 
(Equation  (1.68))  and  at  the  critical  points  i  =  0.  As 
was  previously  mentioned  other  methods  must  be  employed 
to  accurately  determine  the  head  wave  amplitude  near 
this  point. 

When  computing  synthetic  seismograms  in  plane 
layered  media  without  lateral  inhomogeneities  it  is 
convenient  to  use  the  concepts  of  kinematic  and  dynamic 
analogues  (Hron  (1972)).  For  the  rays  shown  in  the 
inserts  of  Figures  10  and  11,  using  the  velocity-depth 
structure  of  Table  2,  it  follows  that  each  of  the  three 
possible  rays  with  four  segments  in  both  layers  1  and  2 
have  the  same  kinematic  properties;  that  is,  they  all 
arrive  at  a  given  epicentral  distance  at  the  same  time. 
Consequently  the  travel  time  need  only  be  computed  once, 
and  as  the  geometrical  spreading  depends  only  on  the 
number  of  ray  segments,  it  too  has  only  to  be  computed 
once  for  each  of  the  three  rays  (a  kinematic  group) . 

Also,  since  two  of  the  three  rays  encounter  the 
same  set  of  reflection  and  transmission  coefficients, 

though  not  in  the  same  order,  they  both  have  identical 
dynamic  properties,  that  is,  their  complex  amplitudes 
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Table  1.1 


% 

anis . 

a 

xi 

a 

zi 

pi 

hi 

a 

X2 

a 

z2 

P2 

0 

1.15 

1.15 

2.02 

2.0 

1.55 

1.41 

2.18 

10 

1.21 

1.10 

2.02 

2.0 

1.55 

1.41 

2.18 

20 

1.26 

1.05 

2.02 

2.0 

1.55 

1.41 

2.18 

Velocities  (a  ,  a  ,  a  ,  a  )  are  in  km/sec,  densities  (p,,p9) 

x  ^  z  ^  x  2  2  2  x  z 

3 

are  in  g/cm  ,  and  the  thicknesses  (h^,!^)  are  in  km.  The  degree 

a  -a 
X  z 

of  anisotropy  is  defined  as  -  x  100%. 

z 

Table  1.2 

Layer  ax  az  p  h 

km/sec  km/sec  gm/cc  km 


1.32 

1.05 

2.02 

2.0 

1.76 

1.41 

2.18 

2.0 

2.20 

1.76 

2.30 

— 

HSPACE 


Figure  1.9  The  effect  on  the  amplitude  of  the  reflected 
wave  (r)  and  the  head  wave  (h)  of  varying 
the  anisotropy  in  layer  1  and  0  percent  (a) 
to  20  percent  (b) .  The  elastic  parameters  of 
the  two-layered  model  are  given  in  Table  1. 
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Figure  1.10  The  logarithmic  amplitude  distance  curve 

for  the  first  dynamic  group  in  the  kinematic 
group  having  four  ray  segments  in  each  of 
layers  1  and  2.  The  three-layered  model 
is  specified  by  Table  2. 
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Figure  1.11  The  logarithmic  amplitude  distance  curve 

for  the  second  dynamic  group  in  the 


kinematic  group  having  four  ray  segments  in 
each  of  layers  1  and  2.  Table  2  contains 
the  parameters  of  this  three-layered  model. 
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are  the  same.  This  being  the  case  amplitude  calculations 
need  only  be  carried  out  once  for  a  dynamic  group. 

Thus  the  three  rays  in  the  kinematic  group  can  be 
divided  into  two  dynamic  groups,  the  second  dynamic  group 
in  this  case  having  only  one  entry.  However,  as  the  ray 
parameter  is  the  same  for  all  three  rays,  the  reflection 
and  transmission  coefficients  common  to  both  groups  of 
dynamic  analogues  are  the  same,  and  a  considerable 
saving  in  computer  time  is  realized. 

1.9  Conclusion 

The  theoretical  development  of  the  dynamic  properties 
for  SH  waves  in  plane  layered  transversely  isotropic  media 
using  asymptotic  ray  theory  has  been  presented  with  an 
emphasis  put  on  having  the  final  formulae  in  a  form  that 
are  readily  programable.  Besides  reflected  waves,  head 
waves  are  investigated  within  the  limited  framework  that 
the  asymptotic  expansion  allows.  The  head  wave  amplitude 
requires  the  evaluation  of  the  first  order  term  in  the 
asymptotic  expansion  while  the  reflected  amplitude 
requires  only  the  zero  order  term. 

A  numerical  discussion  of  a  simple  velocity  depth 
structure  was  given  which  showed  the  effects  of  varying 
anisotropy  on  reflection  and  transmission  coefficients 
and  hence  amplitude  distance  curves.  The  merits  of 
employing  kinematic  and  dynamic  analogues  to  reduce 
computing  costs  when  calculating  synthetic  seismograms 
were  cited  along  with  an  example. 
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CHAPTER  2 

SH  WAVES  IN  LAYERED  TRANSVERSELY  ISOTROPIC 
MEDIA  -  A  WAVE  APPROACH 

2.1  Introduction 

The  problem  of  SH  waves  propagating  in  transversely 
isotropic  homogeneous  plane  layered  media  is  discussed 
through  the  use  of  integral  transforms.  This  procedure 
yields  not  only  the  asymptotic  solution  which  is  also 
attainable  using  asymptotic  ray  series  of  geometrical 
optics  (Karal  and  Keller  (1959)  ,  Cerveny  (1972)  ,  Cerveny 
and  Psencik  (1972) )  but  also  allows  for  the  investigation 
of  the  interference  of  the  reflected  and  head  waves  in 
the  vicinity  of  the  critical  point  (point  of  critical 
refraction) .  It  is  in  this  region  that  asymptotic  ray 
theory  breaks  down  or  in  the  least  introduces  consider¬ 
able  error  in  the  displacement  amplitudes. 

Although  exact  solutions  can  be  obtained  using 
integral  transforms,  they  become  unsuitable  for  numerical 
calculations,  and  consequently  the  exact  solutions  are 
approximated  using  the  method  of  steepest  descents. 

Direct  waves  will  not  be  considered  in  the  following  work. 

The  dynamic  properties  of  reflected  and  head  waves 
at  and  around  the  critical  point  have  been  solved  for 
homogeneous  isotropic  media  in  the  case  of  the  coupled 
P-SV  problem  (Cerveny  (1962),  Cerveny  (1965),  Cerveny 


* 
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(1967),  Cerveny  and  Ravindra  (1971)).  The  notation  in 
this  paper  will  be  similar  to  the  above  mentioned  works 
as  will  be  the  contour  of  integration  used. 


2.2  Mathematical  Preliminaries 

Consider  two  perfectly  elastic  transversely  iso¬ 
tropic  homogeneous  media,  j  and  j+1,  in  welded  contact, 
with  a  point  source  of  SH  waves  situated  at  a  point  zq 
above  the  interface  between  the  two  media.  The  interface 
will  be  defined  as  the  z  =  0  plane.  As  SH  wavefronts  in 
transversely  isotropic  media  are  ellipsoids  of  revolution 
it  will  be  further  assumed  that  the  axes  of  anisotropy 
are  aligned  perpendicular  to  the  z  =  0  plane  (Figure  2.1) 
Introducing  cylindrical  coordinates  (r,<j),z),  the  equation 
of  motion  in  the  j-th  (upper)  medium  can  be  written 


(2.1)  A. 

J 


1  _9_ 
r  9r 


r  9u 


A 


9r 


+  B. 
J 


2 

u 


?-> 
9  u 


9  t‘ 


A  .L  (t)  6  (r)  6  (  z— z  o ) 
r 


where  L(t)  is  the  time  dependence  of  the  source,  and  in 

usual  notation  A.  =  A.^  and  B.  =  A (Cerveny  and 

j  44  j  66 

Psensik  (1972),  Daley  and  Hron  (1977)).  Both  quantities 
Aj  and  B ^  have  dimensions  of  velocity  squared  and  in 
most  seismological  applications  the  condition  A_.>Bj 
holds.  The  plane  of  incidence  can  be  chosen  without 
loss  of  generality  as  the  (j>  =  0  plane  and  thus  the  dis- 

— >■  A  A 

placement  vector  can  be  expressed  as  u ^  =  u^i^  where  i^ 
is  a  unit  vector  in  the  (J)  direction  or  in  the  case  of 


‘ 


Figure 


2.1 


Geometry  of  incidence  and  notation  used  in 

the  case  of  a  point  source  of  SH  waves  at 

r=0,  z=z  .  At  the  interface  at  z=0  the  media 
o 

j  and  j+1  are  in  welded  contact. 
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<J>  =  0 ,  the  y  direction. 

Snell's  Law  is  valid  at  the  interface,  but  if 
expressed  in  terms  of  ray  angles  and  velocities  as 
opposed  to  wavefront  normal  angles  and  velocities,  takes 
on  a  slightly  different  form  (Sato  and  Lapwood  (1968), 
Daley  and  Hron  (1978)).  In  general  the  wavefront  normal 
angles  and  velocities  are  different  from  ray  angles  and 
velocities  (Musgrave  (1970) ) .  The  problem  under  con¬ 
sideration  here  can  be  treated  totally  using  only  ray 
angles  and  velocities;  these  being  the  most  practical, 
for  it  is  along  the  ray  that  the  energy  is  propagated. 
The  statement  of  Snell's  Law  is 


(2.2) 


a  .  sina .  a .  , , sina .  , , 
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ax  .  =  (A .  )  2  - 
3  3 


the  velocity  along 
of  the  ellipsoid, 
the  velocity  along 
of  the  ellipsoid. 


the  major  (r)  axis 


the  minor  (z)  axis 


sin^a . 

=  -i^  + 

ax  . 

3 

_  2 

A .  ax  . 

=  _J_  =  _ 1 

B .  2 

j  az  . 


2 

cos 


2 

az  . 
3 


a . 

1 


defines  the  ray  velocity  a. 


a.  -  the  angle  the  ray  makes  with  the  z  axis. 


- 
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•  _  —  i  CO  t 

A  time  dependence  of  e  will  be  assumed  so  that 

u  =  u  e  and  equation  (2.1)  becomes 


(2.3)  A.  - 


,  ~  r32u.\  32u.  0  6(r)6(z-z  ) 

.  i  i  — - 3_;+  B. - 1  +  u2u.  =  -A. - 

3  r  3r  \  3r  j  3  a,2  3  3  r 


3  z‘ 


with  vector  notation  being  dropped  as  the  result  of  a 
previous  discussion. 

Using  the  Fourier-Bessel  transform  as 


(2.4)  Uj  (k,h,(jo)  = 


r  J  (kr)dr 
o 


-  ,  .  -ihz, 

Uj  (r,z,u))e  dz 


and  its  inverse  as 


(2.5)  u  .  (r ,  z  , 03) 


2  TT 


k  J  (kr)dk 
o 


Uj  (k,h,(o)  e^^dh 


the  solution  of  equation  (2.3)  can  be  written 


(2.8)  u .  =  ^  J  JQ(kr)kdk  j 


exp [ ih ( z-z  ) ] dh 

2  2 
(hZ  +  AZ) 


where  A?  =  i F . k 2  -  )  . 

3  \  3  B. 


Solving  the  integral  over  h  by  the  method  of  residues 
yields 


(2.7)  Uj  =  “F j  e 


7  J  (kr ) exp [ -A . I z-z  | ] kdk 

-iojt  I  o  1  1  o _ 

l  A3 


' 


which  is  the  displacement  for  the  direct  wave  at  any  point 
(r,z)  in  the  j-th  medium. 

The  expressions  for  the  displacement  of  the  reflected 
wave  in  the  j-th  medium  and  the  transmitted  wave  in  the  j+1 
medium  at  an  arbitrary  point  (r,z)  are  respectively 


(2.8)  u  = 


-F  .  e 
3 


-icot 


R.(k)J  (kr ) exp [-A . (z-z  ) ] kdk 
j  o  ^  j  o 

A  . 

3 


and 


(2.9)  u  =  -F . e 
^  D 


-imt 


T.(k)J  (kr ) exp [ A .  , z-A  .  z  ] kdk 
n  o  ^  i  +  l  1  o 


A  . 
3 


where 


A  =  /  f  v  ^ 

j+1  i  j+1  B 


2  '  + 


j+1 


R. (k)  and  T. (k)  are  the  reflection  and  transmission  coeffi- 
:  3 

cients  and  can  be  determined  from  the  interface  conditions 
of  continuity  of  displacement  and  shear  stress.  For 
completeness  these  coefficients  are  discussed  in  Appendix  A, 


2.3  The  Contour  of  Integration 

It  is  convenient  to  express  the  Bessel  function 
JQ(kr)  in  the  equation  for  the  reflected  displacement  in 
terms  of  Hankel  functions  of  the  first  and  second  kind 


as  (Watson  (1944) ) 


- 


- 
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(2.10)  JQ(kr)  =  ± 


H^1}  (kr)  +  H^2)  (kr)  . 


With  the  added  facts  that  H  ^  2  ^  (e  ^kr)  =  -H ^  (kr)  and 

o  o 

Rj (k)  =  Rj (-k)  as  shown  in  Appendix  B ,  it  is  possible  to 
write  the  equation  (2.8)  as 


(2.11)  ur  = 


F  . 

3  -loot 


00  (  1  ) 

r  R. (k)Hv  '  (kr)exp[-A. (z+z  ) ] kdk 


A 


(j0 


With  the  change  of  variable  k  =  -  =  k.q  (2.11)  can  be 


written 


/B .  D 
3 


(2.12)  u  = 
r 


k.F. 
_  3  3 

2 


.  ITT 
-10)t  +  ~2~ 


°r  R.  (q)H^1^  (k  .  qr )  exp  [ik  .  ( 1-F  .  q2  )  ^  (  z+zq)  ]  qdq 

- al - — _ xL _ —  _  - 

'  (1-F.q2)1^ 

OO  D 


On  the  Riemann  sheet  where  the  contour  of  integration 
lies  the  radicals  Aj  and  must  have  the  following 

choices  of  sign  to  ensure  the  vanishing  of  the  exponential 
terms  as  z-*00 


(2.13)  A.  =  -ik. (1-F.q  )  2  =  -ik.w. 

3  3  3  3  3 


Vi  =  _i 


0) 


A 


j  +  1 


1  B~  q  )  lkj+lWj  +  l 


' 
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and  for 


Before  proceeding  further  it  is  useful  to  generalize 
the  problem  of  the  reflected  wave  to  one  of  s-1  trans¬ 
versely  isotropic  layers  overlying  a  halfspace  (the  s-th 
layer) .  The  source  and  receiver  will  be  assumed  to  be  in 
medium  1  just  below  the  first  interface  (the  free  inter¬ 
face)  and  the  distance  between  the  two  is  r.  It  will  be 
assumed  that  there  are  2m  ray  segments  in  the  (s-l)-th 
layer,  implying  that  the  ray  will  be  reflected  m  times 
from  the  s-th  interface.  A  ray  segment  will  be  defined 
as  the  part  of  the  total  ray  between  one  reflection  or 
refraction  and  the  next  (Figure  2.2).  The  modified 
expression  for  the  reflected  displacement  can  then  be 
written  as 


(2.14)  u 


F1k1  — iwt  + 


1 

2 


r 


CO 


m 


(i) 


where 


/I  if  s  =  2 


Z  (q) 


the  product  of  all  reflection  and  trans¬ 
mission  coefficients  encountered  by  the  ray 


. 


, 

. 


Figure  2.2 


A  ray  diagram  showing  the  labelling  of  the 


layers  and  interfaces. 


CN  CO 


CN  _ 

I  I 

to  to  co 


-/ 
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Rs(q) 


except  for  those  m  reflection  coefficients 
at  the  s-th  interface  (s>2) . 
the  reflection  coefficient  at  the  s-th 
interface 


F(q)  =  X  Njhj  (ir, 

3=1  J  \  3, 


B,  \h  t  A. 


1  -  J.  q2 


B 


1 


Nj  -  the  number  of  ray  segments  in  the  j-th  layer 


hj  -  the  thickness  of  the  j-th  layer. 


B. 


In  writing  (2.14)  it  has  been  assumed  that  —  is 

A1 

B1 

the  maximum  value  in  the  sequences  of  values  - —  (j  =  l,s)  . 

A1 

If  this  is  not  the  case  the  radical  in  the  denominator 

l  \  i\h 

in  (2.14)  must  be  replaced  by  I  1  -  q  j  where  k 

1bi  1 

refers  to  the  layer  that  maximizes  —  . 

. 


The  integrand  of  (2.14)  has  2s  branch  points  lying 
on  the  real  axis  of  the  q  plane,  at  the  points 

/Bi\h 

n.  =  ±;  —  i  (j  =  l,s).  The  velocity  structure  will  be 

3  V  j/ 

assumed  such  that  the  minimum  value  of  n.  (j  =  l,s)  is  n  . 

3  s 

This  is  necessary  to  ensure  that  the  first  head  wave  to 
appear  will  be  the  one  from  critical  refraction  at  the 
s-th  interface. 

As  the  integration  path  in  (2.14)  passes  through 
the  point  q=0  where  the  value  of  (k^qr)  becomes 

infinite,  it  will  be  necessary  to  select  a  transformation 


that  will  take  the  contour  from  -°°  to 


oo 


along  the 
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real  q  axis  into  a  new  contour  D  which  does  not  pass 
through  the  point  q=0.  A  discussion  of  this  contour 
will  be  made  shortly. 

If  along  this  new  contour  D ,  k^rlq  |>>1  ,  where 
|q  |  is  the  minimum  value  q  attains  along  D ,  the 
asymptotic  expansion  for  large  argument  of  the  Hankel 
function  (Watson  (1944))  is  valid. 


(2.15) 


H^1*  (kjqrH 


ITT. 


irk^rq 


exp  [ik  rq - £-]  . 


With  (2.15)  equation  (2.14)  becomes 


-icot  + 


177 


(2.16)  u  =  -F1  e 
r  1 


V3* 


27Tr 


TD  -x' 

f  Z(q)Rg(q) exp [ ik1F (q) ] q  2dq 


D 


(1-F1q  ) 


27¥ 


where 


F  (q) 


s=l 

rq  +  \  N . h . 


j  =  l 


D  3 


The  most  common  technique  in  approximating  integrals 
of  the  type  in  (2.16)  is  by  the  method  of  steepest  des¬ 
cents  or  saddle  point  method  because  of  its  simplicity 
and  reasonable  accuracy.  The  saddle  point  is  given  as 
the  solution  of 

Id  F(q)\ 

\  dq  ;q=qo 


(2.17) 


0 


' 
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which  after  some  algebra  yields 


(2.18)  q  = 


az,a1  sman  az.a.sma. 

1  1 _ 1_  1  j 


(j=l,s) 


ax. 


az  . 
1 


If  the  angle  a  is  the  acute  angle  between  the  ray  and  z 
axis  in  the  layer  having  the  second  minimum  value  in  the 

f  B1 

sequence  of  values  —  ,  j  =  l,s  I  ,  then  the  range  of 

\  i 

sina^  is  0<sina^<l  and  consequently  the  range  of  qQ  is 


/  B.  N  :u 


az 


0  ^q 


1 


ax 


Thus  qo  is  always  real  and  positive 


K  j  K 

The  path  of  steepest  descents  as  it  passes  through  the 

3  TT 

saddle  point  makes  an  angle  of  -j-  with  the  positive 

imaginary  q  axis  (Figure  2.3). 

A  basic  assumption  of  the  method  of  steepest  descents 

is  that  the  absolute  value  of  the  exponential  terms  in 

the  integral  decay  rapidly  on  either  side  of  the  saddle 

point  so  that  they  can  be  approximated  by  the  first  term 

in  a  Taylor  series  about  the  saddle  point,  i.e.  their 

value  at  the  saddle  point.  This  approach  is  not  valid 

however,  if  the  non-exponential  terms  vary  rapidly  near  the 

saddle  point.  This  occurs  if  the  saddle  point  is  near  one 

of  the  branch  points,  say  n  ,  due  to  the  radical  wr. 

s  s 

In  the  following  sections,  two  regions  of  the  value 
of  the  saddle  point  will  be  considered: 

(a)  the  saddle  point  qQ  lies  between  0  and  the  first 

branch  point  ng  which  from  now  on  will  be  referred 


- 


Figure  2.3  The  contours  of  integration  for  (a)  the  saddle 
point  lying  between  0  and  the  first  branch 
point  and  (b)  the  saddle  point  lying  between 
the  first  and  second  branch  points.  The  choice 
of  parameterization  forces  the  contour  in  both 
cases  to  pass  to  the  right  of  the  branch  point 
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to  as  n.  The  integral  is  transformed  to  the  new 
contour  shown  in  Figure  2.3  by  the  change  of  variable 

(2.19)  (l-F^q2)55  =  (l-F1qo2)1"2  +  p  exp  f- 

where  p  is  real  (-°o<p<oo)  .  The  integration  of  q  from 
-oo  to  °°  has  been  transformed  to  the  integral  around 
the  partial  semicircle  D  plus  the  integral  around 
the  loop  Do  through  the  saddle  point  qQ  defined  by 
the  parametric  equation  (2.19).  As  the  radius  of 
the  semicircle  is  extended  to  infinity  the  contri¬ 
bution  from  D  goes  to  zero  and  we  are  left  with  the 
integral  along  Dq  which  yields  the  displacement  of 
the  reflected  wave  from  the  s-th  interface. 

(b)  the  saddle  point  qo  is  greater  than  n  but  less  than 
the  value  of  the  next  branch  point  as  in  Figure  2.3 
The  integration  path  now  consists  of  three  parts; 
the  integration  along  the  semicircle  D  which  goes 
to  zero  as  the  radius  tends  to  infinity,  the 
integration  along  Dq  which  gives  the  contribution 
of  the  reflected  wave,  and  the  integration  along 
D*  which  yields  the  contribution  due  to  the  head  wave 
from  the  s-th  interface.  Integration  around  the 
branch  cut  due  to  the  branch  point  n  is  necessary  to 
ensure  that  we  remain  on  the  same  Riemann  sheet. 

The  branch  cut  is  defined  by  the  parametric  equation 


- 


- 
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(2.20)  (1-F1q2)  2  =  (1-F^n2)  2  +  p  exp  (-  real ,  ( 0<p<°°)  . 

Analogous  integration  contours  for  this  problem  for 
isotropic  media  have  been  dealt  with  at  great  length  in 
the  literature  (Cerveny  (1965) ,  Cerveny  and  Ravindra 
(1971),  Ewing  et  al.  (1957),  Brekhovskikh  (1960),  Honda 
and  Nakamura  (1954))  and  for  this  reason  it  would  be 
redundant  to  repeat  the  discussion  here. 

2.4  The  Reflected  Wave 

In  this  section  an  expression  for  the  reflected  dis¬ 
placement,  the  integral  along  contour  Dq,  will  be  obtained 
employing  the  method  of  steepest  descents.  The  exponential 
term  can  be  expanded  in  a  Taylor  series  in  terms  of  the  new 
integration  variable  p  about  p  =  0  (q=q  )  as 

,  ~  2 
k  r  p 

(2.21)  ik^F (q)  =  ik1F(q  )  -  - ^ - 3  + 

2F,  q 
1  ^o 

where 

s-1  B 

(2.22)  r  =  y  H.h.  z—  tana.,  a.  being  the  acute  angle 

]  ]  Bj  11 

the  ray  makes  with  the  z  axis  in  the  j-th  layer. 


With  the  aid  of  equations  (2.2)  and  (2.18)  it  can  be 
shown  that  ik1F(qQ)  =  iooxr  where  xr  is  the  travel  time 
of  the  reflected  wave. 

Substituting  equation  (2.21)  into  (2.16)  the 


, 


‘ 


, 


I 
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expression  for  the  displacement  uq  of  the  reflected  wave 
is  given  by 


(2.23)  u 

o 


-ia)  ( t-x 

r 

e 


Z (p) R™ (p) exp [ -H (q  ) p2] 

b  U 


where 


H<qo> 


Assuming  R  (p)  and  Z (p)  vary  only  slightly  about  p=0 
(q=q  )  i.e.  q  does  not  lie  near  a  branch  point,  the  non¬ 
exponential  terms  in  (2.23)  can  be  approximated  by  their 
value  at  p  =  0  (q=q  ) ,  because  of  the  fact  that  the 
exponential  term  decays  rapidly  in  a  small  region  about 

p  =  0.  It  should  be  remembered  that  as  H(qo)  is  propor- 

/ 

tional  to  frequency  k,  =  - 

V 

the  more  rapid  is  the  decay  about  p  =  0.  Thus  (2.23) 
becomes 


,  the  higher  the  frequency 


(2.24)  u  =  e 


-irn (t-xr)  j  \k 


m 


-  ,  Z  (q  )  R  (q  ) 

2Tirq  I  ^o  s  ^o 


exp [ -H (qQ) p  ] dp 


-io>(t-T  ) 

e  Z(qo)Rs(qo) 


where 


az 


L  = 


1  /  -.rx  \  2 

—  —  is  called  the  geometrical  spreading 


a^  sina^ 


of  the  ray.  Equation  (2.24)  is  the  high  frequency  or 


■V. 


' 

. 
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asymptotic  form  of  the  reflected  wave  which  is  also 

attainable  using  asymptotic  ray  theory. 

In  the  case  where  qQ  lies  close  to  a  branch  point  n, 

approximating  R  (p)  by  R  (q  )  is  not  valid  since  R  (p) 

s  s  o  s 


To  accommodate  for 


this  the  integral  in  (2.23)  is  divided  into  two  parts  by 

expressing  R  (p)  in  the  form 
s 

Rs (P)  =  cx (P)  "  C, (P)ws. 

Explicit  expressions  for  C^(p)  and  (p)  are  given  in 
Appendix  B.  It  follows  that 


R 


s 


For  q  close  to  n,  w  is  approximately  zero  so  that 


o 


s 


C2Cp) 


<<1.  If  this  condition  is  satisfied  it  can  be 


justified  to  approximate  R  m(p)  by  the  first  two  terms 

s 

in  a  binomial  expansion,  viz. 


(2.25)  Rsm(p)  *  C-^Cp)  -  m  Cim  1(p)C2(p)ws 


so  that  equation  (2.23)  can  be  written  in  the  form 


, 
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(2.26)  u  = 


iw  (t-xr)  I  k-^f 


27Tr  i 


v- 


Z  (p)C  m(p)exp[-H(q  )p2]  ^ 
1  °  q2 


J 

•  CO 


Z(p)mCim  1 (p)C2 (p)wgexp [-H (qQ)p2]  ^ 

q 2 


As  the  functions  Z  (p)  ,  C^(p)  and  (p)  do  not  contain 

the  radical  w  ,  they  do  not  change  rapidly  in  the  neigh- 

s 

borhood  of  p  =  0  and  as  an  initial  approximation  it  is 
reasonable  to  assign  to  them  their  value  at  p  =  0.  Thus 


-ioo(t-T  ) 
r 


(2.27)  u  = 
o 


Z(qJ  cim(qo)-mCim  1(qjL 


k^  \  "2 


2urq 


x 


wgexp [-H (qQ) p  ]dp^. 


The  radical  wg  can  be  alternately  written  using 
equation  (2.19)  in  the  form 


/  A  /  2 

(2'28)  «s  =  l1  -  =  t1  -  V 


1 

"X 


F1  2n 


(1-F1n2)^  -  ( 1-F xq 


23g 


2  1 


2 


(1-F^n  )  2+(l-F1q  ) 


r 


| (l-F1nZ)^-(l-F1qo)  2-p  exp 


F,  n  L 


-ITT 


2  1 


2 


(1-F  n  )  2  +  (I-F^q  )  +  P  exp 


-ITT  \  I  % 


J 


X 
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Assuming  that 


2  x  ^  .  ,,  _  2  1" 


(2.29)  H(qQ)2  (1-F1n^)  2  +  (1-F-^  ) 


L. 


>>1 


then  w  can  be  written  approximately  as 
s 


(2.30)  w  =  — i— 
s  „  h 

n 


2  i 


2.35 


( 1-F  n  )  + ( 1-F..  q  ) 
1  1  o 


2.  ^  /-iTr  N~1  3g 

•'l  A  1  '■  • 


( 1-Fn  n  )  - 


-d-F^)2  -p  exp^J 


J 


The  argument  of  the  first  radical  in  the  []  brackets  is 
always  zero  while  that  for  the  second  radical  is,  for  p=0 


2  a- 


2 


(2.31)  arg  :  (1-F^n  )  -(1-Fnq^  )  -p  exp 


r  tt 

/  •  „rii  “  o  for  q  <n 
/  ~ ITT  '  ■  2  nO 


L 


l”o 


\  ~A  f  -  \ 

j  j^O  for  q^>n 


The  values  of  the  argument  are  chosen  to  agree  with 
equation  (2.13)  in  Section  1. 

h 

With  the  change  of  variable  p->H(qQ)  p  and  the  use  of 
equation  (2.30)  the  expression  for  the  reflected  dis¬ 
placement  near  the  branch  point  n  becomes 


(2.32)  u  = 


— ico  ( t— t  )  r 

e  r  „  ,  .  ^  m,  .  im^.m-1,  . 

—  L -  <V  -  —ST  C1  CV 

L  ri 


k  \3g 

C0  (q  )  L  v  H  (q  ) 
2  Mo  2rq  ^o 

\  °/ 


i  r  .-3/4  XC.  _  2  x 


( l-F1n‘ )  s+  ( 1-F^2 )  / 55 


•> 


gl(po) 


where 


‘ 


. 
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CO 


(2.33)  g1(pQ)  =  7T 


—  CO 


PQ-P  exp 


exp (-p2) dp 


and  from  equation  (2.31)  the  argument  of  the  radical 
under  the  integral  in  (2.33)  is 


for  p  <0  (q  <n) 
*0  Mo 


(2.34) 


0  for  p  >0  (a  >n) 


The  frequency  dependent  quantity  |p  |  can  be  thought  of 
as  a  measure  of  the  distance  from  the  critical  point; 
the  critical  point  being  defined  as  the  receiver  location 
at  the  first  interface  where  the  ray  under  consideration 
arrives  for  qQ  =  n.  At  the  critical  point  pQ  =  0 . 

Multiplying  and  dividing  the  second  term  of 
equation  (2.32)  by  pQ  and  employing  the  previously 


defined  relations  for  L  and  w  we  are  left  with 

s 


-ito  (t-x  ) 


(2.35)  u 


e 


■ 


o 


L 


where 


CO 


(2.36)  y-L(P0 


•* 


' 


Graphs  of  the  functions  (modulus  and  phase)  of  9^(P0) 
and  y^(pQ)  are  shown  in  Figure  2.4. 

Equation  (2.35)  can  be  further  simplified  if  it  is 
recalled  that  for  q  -n 


(2.25)'  Rsm(qo> 


^  m  .  , 

C1  (V 


mC 
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and  as  a  result 
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(2.38)  N  =  -raC1m_J-(qo)C2(qo)  1  -  -\'j  (V;L  (pQ)  -1)  . 


Finally  by  setting  qQ  =  n  and  multiplying  and  divid¬ 
ing  the  expression  for  N  by  pQ  with  the  added  facts  that 

C1(n)  =  1 

r  ~ 

(n)  =  -  =  T  which  will  be  called  the  modified 

2  az 

s 

head  wave  coefficient, 

N  can  be  written 


(2.39)  N 


(l-F^n2)®4 


gl(po> 


=  -g* (n) 


with  the  result  that  near  the  critical  point  the 
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The  modulus  and  argument  (in  radians)  of 
y^  and  g^  vs.  p  . 


Figure  2.4 
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expression  for  the  reflected  wave  has  the  form 


-ia)(t-T  ) 
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(2.40)  u 
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where 


g*  (n) 


2  \ 

mr(l-F  n  )  4 


2 . 5  The  Head  Wave 

For  the  region  of  q  referred  to  in  case  (b)  in 
Section  2  there  are  two  displacement  contributions  to  be 
considered;  the  reflected  one  corresponding  to  the 


integral  along  Dq  and  the  head  wave  displacement  corres¬ 


ponding  to  the  integral  around  the  branch  cut  D*.  As  the 
first  part  was  solved  for  in  Section  3  only  the  contribu¬ 
tion  due  to  the  integral  over  D*  will  be  considered  here. 

Head  waves  can  only  arise  at  an  interface  if  the 
incident  wavefront  has  a  finite  radius  of  curvature  in 
the  plane  of  incidence.  (Plane  wave  incidence  cannot 
produce  head  waves) .  The  resulting  head  wave  is  a  three 
dimensional  conical  wave  and  hence  has  an  infinite  radius 
of  curvature  in  the  plane  of  incidence.  Consequently, 
any  other  occurrences  of  incidence  by  the  ray  on  the 
interface  which  produced  the  original  head  wave  cannot 
produce  other  head  wave  segments  (Cerveny  and  Ravindra 
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(1971) ) . 

The  expression  for  the  displacement  due  to  the  head 
wave  is  given  by 


-icot+ 


ITT 
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(2.41)  u*  =  -F  e 


k1N^  Z  (q)  Rg  "  (q)  exp  [  ik^F  (q)  ]  q  ^dq 


27rr 


D. 


(1-F-^) 


where  all  the  quantities  in  (2.41)  have  been  previously 
defined,  the  contour  of  integration  is  shown  in  Figure  3 
and  the  asymptotic  expansion  of  the  Hankel  function  was 
assumed  valid. 

The  integration  contour  D*  which  circumvents  the 


branch  cut  due  to  the  radical  w 

s 

the  parametric  equation 


is  given  by 


(2.20)  '  (1-F^q2)1^  =  ( 1-F 


n2)h 


+  P 


/-iTT 
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In  terms  of  this  new  variable,  the  radical  w  , 

o 

which  has  a  different  sign  on  each  side  of  the  cut,  can 
be  written 


(2.42)  w 
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P  + 


2  (1-F1n2 


where 


on  the  left  hand  side  of  the  cut 
on  the  right  hand  side  of  the  cut 


and 
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2 

arg  p  +  2(1-F^n  )  exp 


We  will  proceed  with  the  solution  of  equation  (2.41) 
under  the  assumption  that  the  maximum  value  of  the 
integral  occurs  near  p  =  0,  and  expand  the  function  in 
the  exponent  in  a  Taylor  series  about  p  =  0  (q=n)  so  that 


where  t,  is  the  travel  time  of  the  head  wave  and  r  =  £  +  r. 


h 


The  quantity  r  was  defined  in  equation  (2.22)  and  in  the 
present  situation  is  evaluated  at  q  =  n.  The  distance  the 
ray  travels  along  the  s-th  interface  is  defined  by  £.  It 
should  be  noted  that  l  =  r-r*  where  r  is  the  source- 
receiver  offset  and  r*,  the  critical  distance,  is  the 
horizontal  distance  from  the  source  at  which  the  head 
wave  from  the  s-th  interface  first  appears.  The  above 
quantities  are  defined  in  Figure  2.5. 

As  before  R  m(q)  can  be  expanded  near  q  =  n  as 
s 


(2.44)  Rsm(q)  *  Cim(q)  -  mC^"1  (q)  C2  (q)  wg. 


The  two  terms  in  the  expansion  of  R  m(q)  result  in 


the  expression  for  the  head  wave  displacement  being  the 
sum  of  two  integrals  around  the  branch  cut.  As  all 
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quantities  in  the  first  integral  are  the  same  on  both 

sides  of  the  cut,  the  contribution  from  this  integral  is 

zero.  The  argument  of  w  changes  by  tt  from  the  left  to 

s 

the  right  hand  side  of  the  cut  and  consequently  the 
integral  involving  this  term  is  non-zero.  The  integration 
for  the  contributing  term  is  then  twice  the  integral  over 
p  from  zero  to  infinity. 

The  head  wave  displacement  has  the  form 


-iu)(t-T  )  -  ^  /  k  \h 
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where  H . (n)  = 
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2  3  * 
2FX  n 


As  in  the  case  of  the  reflected  wave  it  will  be 
assumed  that  the  major  contribution  of  the  integral  occurs 
near  p  =  0  (q=n)  so  the  following  approximations  will  be 
made 
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so  that 


(2.47)  u*  =  2e 
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With  the  change  of  variable  p->  (H*  (n)  )  ^p  equation  (2.47) 
has  the  form 
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The  value  p*  is  proportional  to  Z  and  as  such  is 
equal  to  zero  at  the  critical  point  and  increases  with 
distance  from  the  critical  point.  As  a  result  p*  is 
always  real  and  positive.  For  large  p*  the  following 
approximation  is  valid  (see  Figure  5) . 


(2.50)  g2(p*) 
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6  3/2  for  large  p* . 
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Utilizing  this  asymptotic  form  for  g2(p*)  and 
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The  modulus  and  argument  (in  radians)  of 
la 2  and  g2  vs.  p*- 


Figure  2.5 


after  some  algebraic  manipulation  equation  (2.48)  can  be 
written  as 
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(2.51)  u*  =  e  h  2  F  2  - ,  1 
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where 


r  =  az  r 
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Equation  (2.51)  is  similar  to  the  expression  obtained  for 
the  head  wave  expression  when  geometrical  ray  theory  is 
used . 

Near  the  critical  point  £  and  hence  p*  is  small  so 

that  the  approximation  (2.51)  is  no  longer  valid.  Let 
3 /2 

y2(P*)  =  P*  g2(P*)  from  which  it  follows  that 


(2.52)  u* 


Fie 


-iw (t-Th) 
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2 

where  for  large  p* ,  y2(p*)  ~  e  •  Because  y2(p*)->-0  as 

JAO  equation  (2.52)  is  in  an  inconvenient  form.  Multiply- 

3 /2 

ing  and  dividing  (2.52)  by  Lp*  (L  being  the  geometrical 
spreading  evaluated  at  q=n)  the  following  final  form  is 
obtained  for  head  wave  displacement  near  the  critical 
point 


-iw (t-xh) 

(2.53)  u*  =  - - j- -  Z(n)g*(n)  2  3//4  g2(p*) 


where 
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mf 

(2.54)  g*  (n)  =  — 

(k^nrp 

The  term  r  in  the  denominator  of  equation  (2.54)  should 
actually  be  r,  but  since  r  =  r  +  £  and  i  is  small  near 
the  critical  point  r  has  been  replaced  by  r. 

At  this  point  it  should  be  mentioned  that  (2.40) 
and  (2.53)  are  approximations  for  the  reflected  wave  and 
head  wave  near  the  critical  point,  and  their  accuracy 
diminishes  with  increasing  distance  from  that  point. 
However,  at  large  distances  from  the  critical  point  the 
asymptotic  expressions  (2.24)  and  (2.51)  are  generally 
good  approximations. 

2.6  Interference  Reflected  Head  Wave 

For  the  range  of  qQ(0<qo<n),  only  the  reflected  wave 
from  the  s-th  interface  exists.  The  point  qo=n  corres¬ 
ponds  to  the  onset  of  the  head  wave.  In  Sections  3  and  4 
expressions  were  obtained  for  the  reflected  wave  and  head 
wave,  equations  (2.40)  and  (2.53)  which  are  valid  in  the 
range  (n<qQ<n+e) ,  where  e  is  small  and  its  magnitude 
dependent  on  frequency.  For  values  of  qQ  greater  than  n+e 
the  asymptotic  forms  for  the  reflected  and  head  waves  can 
be  used.  The  sum  of  the  reflected  and  head  wave  displace¬ 
ment  components  in  the  range  (n<qQ<n+e)  will  be  called 
the  interference  wave. 

From  equations  (2.40)  and  (2.53)  the  interference 
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displacement  component  has  the  form 

-ioo(t-T  ) 

r  , 

( 2  •  55)  u  =  u0  +  u*  =  - - l -  Z  (qo)^Rsm(qo)  “  q*(n) 

[i23//4g1  (pQ)  -  i23//4po'^] +g*  (n)  2  3//4exp  (-ip)  g2  (p*)J> 

where  p  =  w"(Tr-T^)  2  and  it  has  been  assumed  that 
Z(qQ)-Z(n).  After  some  rearrangement  (2.55)  becomes 

-ioj(t-T  ) 
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(2.56)  u  =  - - - - Z (qQ) J Rsm(qQ) -g* (n)  [i23/4g1 (po)  - 

•  03/4  ~-3/4  .  ,  .  .-2.  '[ 

-  12  po  -  2  g2  (p* )  exp  (-ip  )]j. 


As  the  quantities  pQ,  p*  and  p  are  all  measures  of 
the  distance  from  the  critical  point,  and  in  the  vicinity 
of  the  critical  point  are  all  approximately  equal,  they 
can  all  be  replaced  by  one  value.  The  quantity  p  is 
chosen  because  of  its  obvious  physical  significance. 

Let 


(2.57)  G  (p)  =  23//4ig1  (p) -23^4ip  ^-2  3//4g2  (p)  exp  ( -ip  2 ) 
so  that 
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With  the  aid  of  equations  (2.33)  and  (2.49)  the 
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expression 

(2.59)  G (p) 


for  G (p)  can  be  written 
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where  arg  p^  =  it  for  p^>0.  Extending  the  limits  of 
integration  of  the  second  term  in  G(p)  results  in  the 
following 
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where  arg  p2  =  0  for  p>0.  (Cerveny  and  Ravindra  (1970)). 

The  form  of  G(p)  in  equation  (2.60)  is  such  that  it 
can  be  expressed  in  terms  of  the  Weber-Hermite  function 
(parabolic  cylinder  function)  Di(p(i-1))  as 


(2.61)  G(p) 


(p(i-l) )-2 


with  D-l  (p(i-l))  being  defined  by  the  identity 
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Figure  2.6  Real  and  imaginary  parts  of  the  Weber  function 
for  negative  and  positive  arguments. 
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00 


(2.62) 


s  exp  -s 


ds 


The  behaviour  of  the  Weber-Hermite  function  is  shown 
graphically  in  Figure  2.6.  Other  properties  of  the  Weber- 


Hermite  function  and  techniques  for  its  numerical  evalua¬ 


tion  can  be  found  in  the  paper  by  Marks  and  Hron  (1977) 
and  a  numerical  tabulation  of  it  is  given  in  Kireyeva 
and  Karpov  (1961) . 

* 

2.7  Numerical  Results 

In  this  section  the  previously  derived  equations  are 
applied  to  a  simple  hypothetical  model  of  two  plane 
layers  overlying  a  half space.  The  description  of  the 
media  is  given  in  Table  1.  Modulus  of  the  amplitude  vs. 
distance  curves  are  plotted  in  the  vicinity  of  the 
critical  point  for  the  two  rays  shown  in  the  inserts  of 
Figures  2.7  and  2.8;  that  is  for  rays  having  multiplicities 
of  one  and  three.  The  significance  of  the  factor  m  in 
the  expression  for  the  head  wave  amplitude  is  that  for  a 
ray  with  multiplicity  m  there  are  m  rays  arriving  at  the 
receiver  with  the  same  kinematic  and  dynamic  properties. 

As  noted  in  the  figure  captions  the  three  curves 


presented  are  for  the  asymptotic  solutions  of  the 
reflected  wave  (a)  and  the  head  wave  (b)  using  equations 
(2.24)  and  (2.51)  and  the  interference  wave  (c)  from 


TABLE  2.1 


Specifications  of  the  media  used  in  Figures  2.7  and  2.8 


Layer  ax  az  p  thickness 

km. /sec.  km. /sec.  gm./cc  km. 
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Figure  2.7  Log.  amplitude-distance  curves  of  the  ray 
shown  in  the  insert  for  (a)  the  reflected 
wave  (b)  the  head  wave  (asymptotic  form) 
and  (c)  the  interference  ref lected-head  wave. 
The  critical  distance  is  denoted  by  r*  and 
the  interference  zone  for  the  source  pulse 
given  in  the  text  lies  between  the  v's. 

It  should  be  remembered  that  curve  (c)  is 
only  valid  in  the  interference  zone. 
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Figure  2.8  Log.  amplitude-distance  curves  for  a  ray  of 
multiplicity  3.  The  logarithmic  singularity 
is  due  to  a  zero  in  the  reflection  coefficient 
of  the  ray  incident  from  below  at  the  second 


interface . 


DISTANCE  -  KM. 
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from  equation  (2.58). 

In  the  asymptotic  limit,  the  head  wave  amplitude  is 
inversely  proportional  to  (r-r*) ,  so  that  at  r=r*  this 
amplitude  becomes  infinite  and  is  obviously  non-physical, 
Also  as  was  previously  mentioned,  since  the  difference 
in  travel  times  between  the  reflected  and  head  waves  is 
small  (i.e.  the  difference  in  travel  times  is  less  than 
the  length  of  the  source  pulse)  in  the  vicinity  of  the 
critical  point,  the  two  arrivals  interfere  with  one 
another  so  that  a  simple  summation  of  the  two  asymptotic 


forms  ut  =  e 


-ioo  (t-T  ) 


-ico  (t  -t,  )  v 

A  +  e  r  A,  in  this 

,  r  h 

L  _  / 


region  introduces  further  error.  In  the  case  under 


consideration,  for  a  source  pulse  L(t)  =  sin(2iTfo) 

/  /  2 Trf  2  \ 

exp  -  - —  ,  with  f  =  20  Hz  and  y  =  4.0  the 

\  .  Y  /  /  ° 

length  of  the  interference  zone  is  approximately  1.9  km 
and  is  denoted  in  Figures  2.7  and  2.8  as  the  distance  between 
the  two  v's.  The  development  of  the  interference  formula 
is  an  attempt  to  reduce  the  aforementioned  inaccuracies 
in  this  zone. 

2 . 8  Conclusion 

The  mathematical  development  of  wave  propagation  for 
SH  waves  in  a  plane  layered  transversely  isotropic  homo¬ 
geneous  media  using  integral  transform  methods  has  been 


1  . 
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. 


. 


presented.  Asymptotic  forms  of  the  reflected  and  head 
wave  displacement  vectors  have  been  derived  along  with 
special  treatment  of  the  interference  of  these  two  waves 
in  the  vicinity  of  the  critical  point,  something  which 
is  not  attainable  if  asymptotic  ray  theory  is  employed. 


. 
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CHAPTER  3 

RAY-REFLECTIVITY  METHOD  FOR  SH  WAVES  IN 
STACKS  OF  THIN  AND  THICK  LAYERS 

3.1  Introduction 

Several  methods  including  plane  wave,  reflectivity, 
integral  transform  and  ray  approaches  have  been  employed 
in  reflection  studies  of  plane  layered  media  to  produce 
synthetic  sections  to  aid  in  the  interpretation  of 
observed  seismograms.  As  would  be  expected  each  of  these 
methods  has  its  own  peculiar  list  of  advantages  and 
disadvantages . 

The  plane  wave  solution  (for  example,  Wuenschel 
(1960) ,  Sherwood  and  Trorey  (1965) ,  Treitel  and  Robinson 
(1966),  Frasier  (1969))  which  is  a  direct  consequence  of 
the  work  of  Thomson  (1950)  and  Haskell  (1953)  has  been 
used  extensively  and  gives  a  good  indication  of  the 
reflection  response  of  a  layered  medium.  However,  a 
plane  wave  solution  is  a  highly  idealized  approximation 
to  an  explosive  point  source  generating  spherical  waves 
as  it  does  not  include  the  geometrical  spreading  of  the 
wavefront  as  it  propagates  through  a  medium.  It  is  also 
instructive  to  be  able  to  identify  individual  arrivals 
on  a  seismic  trace,  an  operation  which  is  difficult  to 
do  using  this  method. 

The  papers  of  Fuchs  and  Muller  (1971),  Fuchs 
(1971)  (1968)  present  the  theory  for  computing  synthetic 
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seismograms  by  the  reflectivity  method  which  is  extremely 
useful  in  refraction  studies  when  investigating  the 
effects  of  a  single  transition  zone.  However  for  sub- 
critical  (reflection)  work  the  numerical  integration 
used  for  refraction  studies  is  not  necessary  and  a 
stationary  phase  approximation  suffices;  whereupon 
inaccuracies  in  amplitudes  are  introduced  if  some  layers 
in  the  transition  zones  are  thick  or  if  multiple  transi¬ 
tion  zones  separated  by  thick  layers  are  to  be  considered. 

Ray  methods  (Helmberger  and  Morris  (1969) ,  (1970) , 

Muller  (1971),  Hron  and  Kanasewich  (1970))  are  useful  in 
that  they  incorporate  geometrical  spreading  and  enable 
individual  arrivals  to  be  identified.  If  a  partial  ray 
expansion  to  the  total  wave  field  is  desired  (Hron  (1971) ) 
many  rays  must  be  considered  to  include  all  multiples  and 
this  places  a  severe  restriction  on  the  number  of  layers 
which  can  reasonably  be  considered.  Albeit  some  success 
has  been  seen  by  including  only  the  contributions  from 
primaries  and  surface  multiples,  the  seismograms  produced 
may  be  misleading  for  some  velocity  depth  structures. 
Another  drawback  of  ray  theory  is  that  it  breaks  down  if 
the  dimensions  of  the  layers  are  of  the  order  of  the 
wavelength  (Cerveny  (1977) )  . 

The  object  of  this  paper  is  to  present  a  compromise 
among  the  methods  mentioned  above  that  will  provide 
greater  flexibility  in  dealing  with  velocity  depth 
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structures  than  any  one  method  alone  and  at  the  same 
time  maintain  or  increase  the  accuracy  of  the  results. 

The  techniques  employed  here  examine  the  velocity 
depth  structure  and  separate  it  into  thick  layers 
separated  by  thin  layered  transition  zones.  The  total 
wave  field  will  be  approximated  by  a  partial  ray 
expansion  in  the  thick  layers  only  while  the  thin 
layered  (transition)  zones  that  have  been  introduced 
will  be  treated  as  quasi-interfaces  and  analogues  of 
reflection  and  transmission  coefficients  are  calculated. 

The  reflection  and  transmission  effects  of  these  thin 
layered  zones  are  referred  to  in  the  literature  as 
reflectivities  and  transmittivities .  Their  computation 
follows  from  the  Thomson-Haskell  technique  and  contain 
the  same  inherent  numerical  problems  (Knopoff  (1964), 

Dunkin  (1965),  Cerveny  (1974)). 

In  the  text  the  simple  case  of  SH  waves  propagating 
in  the  previously  described  medium  will  be  discussed. 

This  simple  case  was  chosen  because  it  presents  the 
basic  idea  of  the  method  without  involving  excessive 
mathematics  which  would  be  necessary  if  the  coupled 
P-SV  case  were  considered.  However,  the  method  is 
readily  applicable  to  the  P-SV  problem  and  has  been 
treated  in  a  relatively  inaccessible  book  in  Russian 
by  Ratnikova  (1973) . 


* 
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3.2  Theoretical  Preliminaries 

Consider  a  system  composed  of  m  plane  isotropic, 
homogeneous  thick  elastic  layers  between  two  half spaces, 
the  upper  of  which  is  a  vacuum  and  the  lower  one  isotropic 
and  homogeneous  (Figure  1).  A  point  source  F(t)  which 
radiates  spherical  SH  waves  is  located  (in  cylindrical 
coordinate  notation)  at  r  =  0,  z  =  0,  4>  =  0 ,  and  a 
receiver  is  situated  at  G,  a  distance  r  away.  The  coor¬ 
dinate  location  z  =  0  corresponds  to  the  free  interface 
between  the  vacuum  and  layer  1  and  z  is  assumed  positive 
downwards . 

The  displacement  vector  of  an  SH  wave  propagating 
in  an  isotropic  homogeneous  medium  is  perpendicular  to 
the  plane  of  incidence  of  the  wave  and  is  decoupled  from 
the  P-SV  wave  propagation.  If  the  plane  of  incidence 
is  chosen  as  the  <f»  =  0  plane,  the  displacement  vector 
can  be  written  without  loss  of  generality  as 

(3.1)  u(r,z,t)  =  u(r,z,t)n 

4> 

where  n  is  a  unit  vector  in  the  <J>  direction,  and  the 

4> 

Fourier  time  transformed  disturbance  due  to  an  arbitrary 
ray  propagating  in  some  or  all  of  the  m  layers  at  the 
point  G  can  be  written  as  (Daley  and  Hron  (1979)) 


oo 


_ 


. 
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Figure  3.1 


Geometry  of  the  media  for  (a)  plane  inter 
faces  between  thick  layers  and  (b)  quasi¬ 
interfaces  between  the  thick  layers 
composed  of  thin  layered  transition  zones 


' 
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where  L(a))  is  the  Fourier  transform  of  the  source  pulse 
F(t),  J^(kr)  is  the  Bessel  function  of  the  first  kind  and 
zero  order,  R(k)  is  the  product  of  all  the  reflection  and 
transmission  coefficients  encountered  by  the  ray,  and  the 
quantities  are  defined  by 


(3.3) 


v  = 
n 


(k2  -  k2 ) ^ 
n 


■  ,  i,2  .  2  |  . 

-i  k  -  k  )  2 
'  n  1 


c>  8 


n 


c<£ 


n 


The  horizontal  wavenumber  k  is  given  as  k  =  —  and 

c 

k  =  ,  £  is  the  shear  velocity  in  the  n-th  layer, 

n  p  n 

n 

a)  is  circular  frequency  and  the  quantities  and  h^ 
are  the  number  of  ray  segments  in  the  n-th  layer  and 
the  thickness  of  the  n-th  layer  respectively. 

It  will  now  be  assumed  that  some  or  all  of  the 
interfaces  encountered  by  the  ray  are  replaced  by  thin 
transition  zones;  thin  in  the  fact  that  their  thicknesses 
are  negligible  when  compared  to  the  thicknesses  of  the 
other  thick  layers.  These  transition  zones  may  take  the 
form  of  a  single  thin  high  or  low  velocity  channel  or 
stack  of  several  thin  layers  of  varying  velocity  (Figure  lb). 

In  this  case  equation  (3.2)  is  changed  in  that  the 
frequency  independent  term  R(k)  must  be  replaced  by  the 
product  of  the  now  frequency  dependent  reflectivities  and 
transmittivities  P(oo,k)  so  that 


. 


' 


89 


CO 


(3.4)  u  (r ,  0  ,  go)  =  L  (go) 


0 


r  m 

P  (go  , k)  exp  !  -i  T  N  h  v 

n  n 

L  n=l 


n  n  n 


One  of  the  most  common  methods  of  evaluating  integrals  of 
the  type  in  equation  (3.4)  is  by  the  stationary  phase 
approximation.  By  comparison  with  direct  numerical 
integration ,  this  method  will  be  shown  to  be  a  quite 
acceptable  approximation  for  small  source-receiver  offsets. 
The  condition  which  must  be  satisfied  for  the  successful 
implementation  of  the  stationary  phase  approximation  is 
that  the  non-exponential  terms  in  equation  (3.4)  are 
slowly  varying  implying  that  the  region  of  applicability 
must  be  removed  from  all  singular  points  (branch  points 
and  poles) . 

3.3  Reflectivity  and  Transmittivity 

For  completeness,  the  reflectivity  and  transmittivity 
of  a  thin  layered  or  transition  zone  will  be  briefly  dis¬ 
cussed  in  this  section.  Similar  treatments  of  this 
problem  may  be  found  in  the  literature,  especially  Cerveny 
(1974)  and  Molotkov  et  al  (1972). 

A  thin  layered  zone  composed  of  s  isotropic,  homo¬ 
geneous  elastic  layers  will  be  assumed  to  lie  in  welded 
contact  between  the  thick  layers  (layers  0  and  s+1)  which 
are  also  isotropic  homogeneous  and  elastic.  The  vertical 
coordinate  designating  the  bottom  of  the  n-th  layer  will 
be  denoted  as  z  and  the  thickness  of  the  n-th  layer  as 


n 
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hn*  Without  loss  of  generality,  let  zq  =  0.  Each  layer 

can  be  fully  defined  for  SH  propagation  by  the  density 

y 

p  and  the  shear  wave  velocity  3  =  (y  /p  )  2  when  y  is 

n  2  n  n  Kn  n 

the  shear  modulus. 

A  slight  deviation  in  notation  will  be  used  in  this 
section  as  it  is  convenient  to  assume  that  a  plane  wave 
impinges  on  the  thin  layered  zone  when  deriving  the 
expressions  for  the  reflectivity  and  transmittivity . 

This  will  require  the  introduction  of  a  Cartesian  coor¬ 
dinate  system,  with  the  plane  of  incidence  being  defined 
as  the  (x,z)  plane  (Figure  2).  This  assumption  does  not 
affect  the  results  as  it  has  been  shown  by  Fuchs  (1968) , 
(1971)  that  the  kernel  P(m,k)  of  the  integral  in 
equation  (3.4)  is  the  same  whether  a  spherical  or  plane 
wave  approach  is  used. 

a2 

2  1  9  Un 

The  solution  of  the  wave  equation  V  u  =  — ~ 

n  en2  at2 

in  each  of  the  layers  has  the  form 


(3.5) 


u  =  cf>  +4) 
n  n  n 


where  the  "+"  and  signs  indicate  downward  and  upward 

propagating  plane  waves  and  have  the  form 


(3.6) 


exp[  +  i  (z-z  1)  vn]  exp[i  (cot-kx)  ] 


where  as  before 


' 


Figure  3.2  Description  of  layers  and  layer  parameters 


for  computing  the  reflectivities  and 

transmittivities  of  a  thin  layered  zone. 

3  and  p  denote  the  shear  wave  velocity 
n  n 

and  density  respectively  of  the  n-th  layer. 
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r 


(3.7)  vn  = 


2  2  -5 

(k  -  kz)  2 
n 


2  ,  2 


!  — i  (  k  -  k  ) 

l  n 


for  c^8 


n 


for  c<8  . 
n 


The  quantities  A*  and  A^  are  complex  functions  of  w 

and  k  to  be  determined  from  the  boundary  conditions. 

Introducing  the  column  vectors  $  and  S 

n  n 


(3.8)  $  = 

n 


1  c 

-e- 

1 _ 

u 

n 

S  = 
n 

4>+ 

Yn 

(a  ) 
yz  n 

3u 


where  (a  )  =  y  „ 

yz  n  n  dz 


n 


is  the  shear  stress,  the  following 


expression  can  be  derived 
(3.9)  Sn  =  Tn*n 


where 


r 


(3.10)  Tn  = 


iy  v 

n  n 


-iy  v 
n  n 


and 


i  iy  v 
n  n 


(3.11)  T 


-1 

n 


2  iy  v 
Kn  n 


n  n 


1 


-1 
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It  also  follows  from  equation  (3.9)  that 

(3.12)  $  =  T_1S  . 

n  n  n 

The  vector  at  the  lower  boundary  of  the  n-th 
layer  can  be  related  to  at  the  upper  boundary  of  the 
n-th  layer  as 


(3.13)  $  ( z  )  =  E  $  (z  _ ) 

n  n  n  n  n-1 


where  it  can  be  deduced  that 


(3.14)  L  = 
n 


exp  iQ  0 

n 


0 


exp  -iQ 


n- 


with  Q  =  (z  -z  n)v  =  h  v  . 
n  n  n-1  n  n  n 

Utilizing  equations  (3.9) ,  (3.12)  and  (3.13)  the 

following  relation  between  the  vector  at  the  bottom 

of  the  n-th  layer  (zn)  and  the  top  of  the  n-th  layer 

(z  . )  can  be  obtained 
n-1 

(3.15)  S(z)=T$(z)=TE$(z  .  )  =  T  L  T_1  S  (z  .  ) 

n  n  n  n  n  n  n  n  n-1  n  n  n  n  n-1 


As  the  boundaries  between  the  layers  are  in  welded 
contact  there  must  be  continuity  of  displacement  (u)  and 


shear  stress  (a  )  at  each  interface,  so  that  at  the  n-th 

yz 


interface 


. 
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(3.16) 


WV  = 


S  (z  ) 
n  n 


Continuing  this  analysis  allows  for  the  expression 
of  the  vector  $  +^,  at  the  boundary  of  the  lower  half space 
(z=zg)  in  terms  of  at  the  boundary  of  the  upper  half¬ 
space  in  the  form 


(3.17)  .  ( z  )  =  D  $  (z  ) 

s+1  s  o  o 


where  the  matrix  D  is  given  by 


(3.18)  D  =  T  • (T  E  T”1)  • (T  ,E  T'1  )  . . .  (T,E  t'1)  -T 

s+1  s  s  s  s-1  s-1  s-1  111  o 


=  T  ]■  «C  *C  ...  .C.,  -T 
s+1  s  s-1  1  o 


and  it  can  be  shown  that  the  matrix  C  has  the  form 

n 


(3.19)  C  =  T  E  T  1 
v  n  n  n  n 


cosQ 


n 


•y  v  sinQ 
n  n  n 


sinQ 


n  i 


y  v 

n  n 


cosQ 


nj 


The  system  of  equations  from  equation  (3.17)  can  be 


written 


(3.20) 


s+1 


Ls+1 


u 


11  D12 


°21  D22 


r  -  "l 

A 

o 


A+ 

o 


from  which  are  obtained  the  relations 
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(3.21) 


A 

o 


s+1 


D 


12 


D 


11 


^21  A- 
D11  S+1 


+ 


det  D 


D 


11 


Introducing  the  notation  R^,  R12,  R22  and  R21  as  in 

Figure  3  and  using  the  relations  given  in  equation  (3.21), 

the  reflectivity  and  transmittivity  of  a  thin  layered  zone 

can  be  obtained.  For  incidence  from  the  upper  thick 

layer  As+p=0  anc^  f°r  incidence  from  the  lower  thick 

layer  A+=0  so  that 
o 


D 


12 


lll 


D 


11 


(3.22) 


R 


det  D 


12 


D 


11 


D 


R 


21 


22 


D 


11 


R 


21  D 


11 


The  term  r(to,k)  in  equation  (3.4)  is  the  appropriate 
product  of  these  reflectivities  and  transmittivities 
evaluated  at  the  transition  zones  encountered  by  the  ray. 
It  is  easily  shown  that  for  a  single  interface  between 
two  layers  equation  (3.22)  produces  ordinary  reflection 
and  transmission  coefficients  for  the  case  of  plane  waves 
incident  on  a  plane  boundary. 


Figure  3.3 


m 


Definition  of  the  notation  used 
describing  the  reflectivity  and  trans- 
mittivity  of  a  thin  layered  zone. 
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3.4  Stationary  Phase  Approximation 


Equation  (3.4)  gives  the  exact  displacement  due  to 
an  arbitrary  ray  arriving  at  G.  However,  preserving  thi 
exactness  requires  that  numerical  integration  be  used  to 
evaluate  the  integral.  This  method  is  unsuitable 
especially  if  the  formula  is  used  to  construct  synthetic 
seismograms  where  a  great  number  of  rays  are  to  be  con¬ 
sidered,  and  hence  a  great  many  numerical  integrations 
must  be  performed.  Consequently,  the  integral  in 
equation  (3.4)  will  be  approximated  by  the  method  of 
stationary  phase  which  as  previously  mentioned  corres¬ 
ponds  quite  reasonably  to  the  exact  solution  for  small 
source-receiver  offsets. 

For  values  of  the  integration  variable  k  greater 

than  some  value  kQ  where  kD  is  the  minimum  in  the 

series  of  values  (kc  ,  kD  . . .  kQ  ) ,  j  being  the  deepest 

61  e2 

thick  layer  the  ray  traverses  (l<j<m),  (l<£<j),  the 

exponent  in  equation  (3.4)  will  have  a  real  negative 
term  implying  an  exponentially  decaying  integrand.  If 
the  thickness  of  this  layer  £  is  large,  the  integrand 
will  decay  fairly  rapidly  for  values  of  k>k^  and  for 


£ 


this  reason  values  of  k>kD  will  not  be  considered. 


$ 


£ 


Thus  after  replacing  the  upper  limit  of  integration  by 


kD  equation  (3.4)  becomes 
3£ 


' 
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(3.23)  u(r,0,co)  =L(co) 


£ 


P  ( 03 ,  k)  exp 


0 


r  3 

-i  7  N  h  v 
n  n  n 

_  n=l 


^  .  kdk 

J  (kr)  - - 

o  IV, 


Introducing  the  change  of  variable 


(3.24)  k  =  k^q,  the  radicals  become 


2  ^  ,  ,  , ,  2  ^ 


(3.25)  vn  =  k^  l-(3n/31q)“  44  with  v1  =  k^l-q") 


and  the  limits  of  integration  are  given  by  (0,3^/3^). 
Equation  (3.23)  can  be  rewritten  as 


(3.26)  u(r,0,oo)  =  k^L(rn) 


0 


r  j 


P(u),q)exp  -i  T  N  h  v  (q) 
:  n=i  n  n  n  ^ 


x  JQ  (k-j^qr ) 


qdq 


2  - 


i ( 1-q  ) 


If  it  is  assumed  that  the  argument  of  the  Bessel 
function,  k^qr,  is  large,  Jo(k^qr)  can  be  used  in  its 
asymptotic  form,  viz. 


(3.27)  Jo(kxqr) 


/2TTk-^qr 


exp 


nr  n 

i(k,qr  +  J)  +  exp  -i(k,qr  +  J 

J  _  *JJ 


Substituting  equation  (3.27)  into  equation  (3.26)  has 


(3.28)  u(r,0,(u)  =  l  J  L  tw)  e 


/..l 


Bl/Bl 


0 


p (W,q) 
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exp'-i  y  N  h  v  (q) — k_ qr 
L,  n  n  n  ^  l^1 


V.  3=1 


k,  \h 


i-T 


<V^ 


27rr 


L  (w)  e 


-i37T/4 


7#, 

q2dq 

'  (1-q2)* 

P  (u),q) 


0 


I  N 


exp.-ij  l  Nnhnvn(q)  +klqr 

L  n=i 


q2dq 


(1-q  ) 


2 


It  will  be  assumed  that  the  integrals  in  equation 

(3.28)  can  be  approximated  by  the  method  of  stationary 
phase  and  further  that  the  phase  of  the  product  of  the 
reflectivities  and  transmittivities  need  not  be  included 
in  computing  the  stationary  point.  This  second  assump¬ 
tion  is  valid  only  for  epicentral  distances  less  than  the 
critical  distance,  as  after  the  critical  distance  the 
phase  of  P(oj,q)  varies  greatly,  while  before  the  critical 
distance  its  phase  is  approximately  constant  for  changing 
values  of  q.  Defining 

(3.29)  f  (q)  =  ±kn  qr  +  T  N  h  v  (q) 

^  1 ^  L.  n  n  n  ^ 

n=l 

it  follows  that 


(3.30) 


df  (q)  _  _ 

dq  -V  klq 


I 


/ 


3 


n=l 


N  h  1 - 

n  n ,  g 


n 


2 


3 


n 


lf  qo  37  1 


1 


In  the  "+"  case. 
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=  tan0  where  6  is  the  acute  angle  a  ray  segment  in  the 
n  n 

n-th  layer  makes  with  the  vertical  axis.  Thus  it  is 
required  that 


(3.31) 


q  =  sin0n . 

1 


The  case  has  no  stationary  point  in  the  range 

31 

of  definition  of  0<q<  -5—  so  that  the  first  integral 

3£ 

in  equation  (2.28)  may  be  neglected. 

Using  the  standard  formula  for  the  stationary 
phase  approximation  to  an  integral,  the  displacement 
can  be  written  as 


(3.32) 


u  (r ,  0  ,  m)  =  P(u),k) 


L(m) 

COS0 


k^sin0  ^ 


rf ' ' (qQ) 


exp [ -if (qQ) ] 


with 


f(qo)  =  ki 


sin0 


m  $ 

,  r  +  J  N  h  ~ 
1  l  n  n  3 

n=l  n 


COS0 


n 


and 


m 


f"(qQ)=ki  I 


N  h  3 
n  n  n 


n=l  cos^0  31 
n 


As  was  previously  mentioned  the  formula  (3.32)  is  valid 
only  for  small  source  receiver  offsets,  so  that  none  of 
the  angles  the  ray  segments  make  with  the  vertical  axis 


approach  critical  angle  values.  Also,  it  is  assumed 
that  the  thickness  of  a  transition  zone  is  small  when 
compared  with  the  so-called  thick  layers,  as  spherical 
divergence  or  geometrical  spreading  is  only  incorporated 
in  equation  (3.32)  for  these  thick  layers.  This  geo¬ 
metrical  spreading  is  given  by  the  k  //rf ' ' (qQ)  term. 

3.5  Numerical  Results 

A  simple  hypothetical  model  with  two  transition 
zones  and  two  thick  layers  (Figure  ,  Table  1)  is 
used  to  compare  three  methods  of  computing  synthetic 
siesmograms;  (a)  ray  theory , 3.4(b)  numerical  integration 
employing  the  reflectivity  method,  and  (c)  the  station¬ 
ary  phase  approximation  discussed  in  this  paper.  The 
object  of  this  comparison  is  to  demonstrate  the 
efficiency  and  reasonable  accuracy  of  the  ray- 
reflectivity  method. 

The  number  of  total  layers  was  kept  small  at  six 
to  enable  a  partial  ray  expansion  (Hron  (1971))  to 
approximate  the  total  wave  field  by  the  ray  method.  In 
this  method  all  possible  rays  up  to  those  with  a  maximum 
of  eighteen  segments  were  generated  in  the  six  layers 
and  included  in  the  seismogram  provided  their  arrivals 
were  within  the  time  window  specified.  In  the  ray- 
reflectivity  method,  rays  needed  only  be  generated  in 
the  two  thick  layers. 


, 


Figure  3.4 


Velocity-depth  structure  of  the  model. 


VELOCITY  KM.  per  SEC. 


Three  separate  programs  were  used  for  each  of  the 
traces  shown  in  Figure  3.5.  The  ray  method  involved 
convolution  in  the  time  domain  while  the  other  two, 
numerical  integration  and  ray-reflectivity,  employed  an 
FFT  algorithm  (Cooley  and  Tukey  (1965))  when  converting 
from  the  frequency  to  time  domain. 

The  source  pulse  used  was 


(3.33)  F(t)  =  sin27Tf  t  exp  - 

° 


£uf  t \  2  \ 
_ o  \ 

Y 


—  °°<  t<°° 


where  f  is  the  predominant  frequency  of  the  pulse  and 
y  is  a  damping  factor.  The  Fourier  transform  of  (3.33)  i 


(3.34)  0  L(w) 


ITT  Y 
2  7T  f 

o 


exp 


/  2 

smh  -y  >- 
47rf 


o  / 


This  particular  source  pulse  was  chosen  because  of  the 
ease  with  which  it  may  be  handled  in  the  frequency  domain 
as  its  spectrum  is  a  Gaussian  curve  whose  maximum  is 
centered  at  f  .  However,  all  three  of  the  programs  used 
can  incorporate  any  source  pulse  whether  in  analytic  or 
digital  form. 

Figure  3.5  indicates  that  the  ray-reflectivity  method 
gives  a  good  approximation  to  the  exact  (numerical  inte¬ 
gration)  solution  even  though  in  terms  of  CPU  time  it 
was  about  forty  times  faster.  The  times  for  the  ray 
approach  and  the  ray-reflectivity  method  were  about  equal 
but  it  must  be  remembered  that  the  model  chosen  was  to 


‘ 
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Figure  3.5  A  comparison  of  three  methods  of  computing 

synthetic  seismograms  for  the  model  shown  in 
Fig.  4.  The  three  methods  are  (a)  asymptotic 
ray  theory  (b)  numerical  integration  and 
(c)  ray-reflectivity  method.  The  notation 
on  the  time  axis  in  cases  (a)  and  (c) 
identify  the  ray  code  of  the  arrival  and  its 
amplitude.  The  source-receiver  offset  in 
all  cases  is  1.0  Kft. 
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facilitate  the  ray  approach  and  that  the  ray-reflectivity 
may  handle  much  more  complicated  structures. 

As  the  ratio  of  the  dimension  of  the  stack  of  thin 
layers  to  the  dimension  of  the  thick  layers  used  here 
is  quite  large  some  error  is  introduced  as  can  be  seen 
when  the  ray-reflectivity  trace  is  compared  to  the  trace 
obtained  by  numerical  integration.  This  error  is  reduced 
as  the  aforementioned  ratio  is  decreased. 

3 . 6  Conclusion 

The  technique  presented  is  well  suited  for  the 
computation  of  synthetic  seismograms  in  oil  exploration 
for  comparison  and  interpretation  of  real  data,  as  it 
correctly  produces  the  extensive  interference  phenomena 
due  to  the  fine  layering.  This  is  done  without  any 
effort  to  decompose  them  into  individual  ray  paths  in 
the  thin  layers.  Physically  justifiable  amplitudes  are 
produced  as  the  method  accounts  for  the  divergence  (in 
the  thick  layers)  of  wavefronts  radiated  from  a  point 
source  and  seismic  energy  partitioning  due  to  reflec¬ 
tions  and  transmission  through  the  thin  layer  stacks. 

The  main  advantage  of  this  method  is  that  a  small 
number  of  rays  may  be  used  to  generate  synthetic  seis¬ 
mograms  in  a  many  layered  structure  and  their  arrivals 
can  be  identified  in  the  seismograms. 
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CHAPTER  4 

A  COMPARISON  OF  SYNTHETIC  SEISMOGRAMS  FOR  A  THINLY 
STRATIFIED  MEDIUM  AND  A  TRANSVERSELY  ISOTROPIC  MEDIUM 

4.1  Introduction 

That  a  periodic  isotropic  two  layered  (PITL)  medium 
behaves  as  a  homogeneous  transversely  isotropic  (HTI) 
medium  in  its  kinematic  and  dynamic  properties  has  been 
shown  to  be  a  valid  assumption  as  far  as  the  propagation 
of  elastic  waves  is  concerned  provided  that  the  seismic 
wavelengths  used  are  large  when  compared  to  the  thicknesses 
of  the  individual  isotropic  layers. 

A  review  of  the  literature  involving  this  problem  to 
1962  is  contained  in  the  paper  of  Backus  (1962).  Notable 
among  the  references  are  the  works  of  Postma  (1955) ,  Uhrig 
and  Van  Melle  (1955) ,  White  and  Agona  (1955) ,  Rytov  (1956) , 
Riznickenko  (1949)  and  Anderson  (1961).  Subsequent  to 
these  are  the  publications  of  Gassmann  (1966)  in  which  is 
presented  a  comprehensive  tutorial  on  the  many  aspects  of 
elastic  wave  propagation  in  (HTI)  media  and  Levin  (1979) 
where  instructional  figures  of  wavefront  profiles  of  (HTI) 
equivalents  of  actual  (PITL)  media  may  be  found. 

However,  to  the  author's  knowledge,  no  graphic 
comparison  of  methods  can  be  found  in  the  literature 
that  these  two  formulations  are  equivalent.  The  intent 
of  this  paper  is  to  consider  a  very  simple  example  of  a 
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(PITL)  medium  composed  of  2N  layers  sandwiched  between  an 
isotropic  layer  and  an  isotropic  halfspace  and  to  compare 
kinematic  and  dynamic  properties  via  synthetic  seismograms 
to  a  similar  geometry  in  which  the  (PITL)  layer  has  been 
replaced  by  its  equivalent  (HTI )  medium  (Figure  4.1). 

This  comparison  will  be  done  only  for  SH  waves. 

The  two  methods  which  will  be  employed  to  compute 
synthetic  seismograms  are  (a)  asymptotic  ray  theory  for 
which  the  required  theory  and  formulae  were  developed 
for  SH  waves  in  an  (HTI)  medium  in  Daley  and  Hron  (1979) 
and  (b)  the  numerical  integration  reflectivity  method  in 
which  each  of  the  2N  layers  of  the  (PITL)  medium  are 
individually  considered. 

4.2  Notation  and  the  Two  Methods 

Aside  from  the  density  p,  two  elastic  parameters 

are  required  to  fully  describe  SH  wave  propagation  in  an 

(HTI)  medium,  in  which  the  wavefront  is  an  ellipsoid  of 

revolution.  The  parameters  are  and  are  such 

C  C 

^44  66 

that  A..  =  -  and  Arr  =  -  have  the  dimensions  of 

44  p  66  p 

velocity  squared.  The  square  root  of  A^  and  A^  are 
respectively,  the  velocity  along  the  horizontal  (r)  and 
vertical  (z)  axis  of  the  ellipsoid  describing  the  wave 
surface. 

The  quantities  and  Cg6  and  density  p  are 

obtained  from  the  shear  moduli  \i1  and  and  densities 

p  and  p  of  the  two  isotropic  layers  comprising  the 


■ 
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(PITL)  medium  by  the  following  averaging  techniques 


(4.1) 


'44 


(dI+dII)yIyII 

dlylI+dIIyI 


(4.2) 


66 


yIdI  +  ^IIdII 

dI  +  dII 


(4.3) 


P  = 


dIpI  +  dIIpII 

dI  +  dII 


where  and  d^  are  the  thicknesses  of  isotropic  layers 
I  and  II  and  N  (d^td^ )  is  the  thickness  of  the  (PITL) 
medium. 

In  what  follows  the  subscripts  1,  2  and  3  will  refer 
respectively  to  quantities  in  layers  1  and  2  and  the 
half space.  In  particular  is  the  shear  wave  velocity 
and  is  the  density  in  the  i-th  layer.  In  the  (HTI ) 
equivalent  of  the  (PITL)  medium  (layer  2)  the  shear  wave 
velocity  (the  velocity  of  energy  propagation)  is 
defined  by 


(4.4) 


1  _  sin  (j)  cos  c{) 

d  2  A..  A  _  ' 

$2  44  66 


4)  being  the  acute  angle  between  the  ray  and  the  vertical 
(z)  axis  (Daley  and  Hron  (1979)). 

As  was  previously  mentioned,  the  asymptotic  ray 
theory  solution  to  the  propagation  of  SH  waves  in  a  (HTI) 


. 


Figure  4.1  Geometry  of  the  media  and  definition  of 


parameters . 
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medium  has  been  presented  by  Daley  and  Hron  (1979)  and  all 
relevant  formulae  for  the  computation  of  synthetic  seismo¬ 
grams  may  be  found  there.  In  that  work  the  first  order 
approximation  of  asymptotic  ray  theory  is  employed  to 
determine  the  solution  for  the  disturbance  at  a  given 
source-receiver  offset  and  is  assumed  to  be  the  sum  of 
the  disturbances  of  a  finite  number  of  rays  which  is 
taken  to  approximate  the  total  wave  field.  As  the 
numerical  integration  method  will  yield  head  wave 
arrivals  as  well  as  reflected  arrivals  the  high  frequency 
ray  approximation  for  head  wave  arrivals  is  also  included 
in  the  seismograms. 

The  time  transformed  expression  for  the  SH  displace¬ 


ment  (u)  at  a  point  (r,z)  in  cylindrical  polar  notation 
within  an  isotropic  layer  due  to  the  reflection  from  a 
stack  of  layers  is 


00 


(4.5) 


R(k,m)  JQ(kr) 


0 


exp  [  iVj(  z-2h)  ] 


kdk  -t 

— —  J 


where 


circular  frequency 


k  -  horizontal  wave  number 


r  -  source  receiver  offset 


h  -  thickness  of  overlaying  isotropic  layer 


' 

, 

■ 


Ill 


J 

o 

L(u>) 

R(k#oj) 


zero  order  Bessel  function  of  the  first 
kind 

time  transform  of  the  source  function  f (t) 
reflectivity  of  the  thin  layered  zone 


Ck2  -  k2)^ 

for 

c>3 

-i(|k2  -  k2\)h 

for 

c<  3 

U) 

c  =  k 

-  unit  vector  perpendicular  to  the  plane  of 
incidence  ((r,z)  plane)  and  the  orientation 
of  axes  and  geometry  is  shown  in  Figure 
(4.1).  If  the  receiver  is  located  on  the 
surface,  z=0. 


The  term  R(aj,k)  is  the  reflectivity  of  the  PITL  layer 
and  is  defined  as 


(4.6) 


R(w,k) 


D(l,2) 

D  ( 1 , 1 ) 


where 


(4.7)  D  =  T3X  (CICII)n  Tx 

and 


(4.8)  Tx  = 

!  . 

L  iy1v1  _1^ivi 
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iy3V3 


(4.9)  T 


-1  =  1 

3  2iy3v3  J 


iy3V3 


-1 


sinQ 


cosQj 


u 


(4.10) 


C.  = 
3 


y  .  v  . 
3  3 


j  =  I, II 


-y^VjSinQj  cosQ^ 


(4.11)  with  v  = 


K 


2  ^ 


-± < i k  -k  i) 


for  c>3 


for  c<3 


K 


K 


k  =  1,3,1,11. 


and  0^  =  ^jvj  '  kein<3  the  thickness  of  the  j-th 

isotropic  layer  making  up  the  PITL  medium. 

To  simplify  computation  and  to  greatly  reduce 
computer  time,  the  matrix  A  =  can  be  diagonalized 

which  minimizes  the  number  of  operations  to  compute  An. 
Thus  A  takes  the  form 

(4.12)  A  =  P-1  A  P 

where  A  is  a  two  by  two  diagonal  matrix  whose  elements 
are  the  eigenvalues  (A^  and  \^)  of  A.  These  eigenvalues 
are  obtained  from  the  solution  of 


. 


113 


*N 


(4.13)  >  -  2cosQ  cosQ__  +  6  sinQ  sinQ  >  X 

^  III  I  1 1  i 

2  2  2  2 
+  Scos  QjCos  QXI  +  sin  Q];sin 

-SsinQjSinQjjCOsQ^cosQjjj  =  0 


with 


.  _  uivi  ,  Miivn 
piivu  yivi 


P  is  the  matrix  whose  columns  are  the  eigenvectors 
corresponding  to  the  two  eigenvalues  and  P  ^  is  its 
inverse.  As  a  consequence 

(4.14)  An  =  P-1AnP. 

The  determination  of  the  eigenvalues  and  ultimately  A 
and  P  is  done  numerically  utilizing  the  IMSL  routing 
EIGCC . 


4.3  Numerical  Discussion 

As  was  indicated  earlier  the  model  to  be  used  in 
comparing  the  two  techniques  in  computing  synthetic 
seismograms  is  an  isotropic  layer  over  a  (PITL)  or 
(HTI )  layer  over  a  half  space.  The  (PITL)  medium  was 
chosen  to  be  composed  of  500  combinations  of  alternating 
layers  of  shale  and  limestone  each  of  one  meter  thickness 
for  a  total  layer  thickness  of  one  kilometer.  The  velo¬ 
cities  and  densities  used  are  given  in  Table  (4.1). 
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In  the  numerical  integration  approach  a  trans¬ 
formation  was  made  so  that  the  integration  variable  was 
the  angle  of  incidence  in  the  overlaying  isotropic  layer 
(Fuchs  1971).  Angular  increments  of  0.1  degree  were 
chosen  and  a  trapezoidal  numerical  integration  scheme 
was  employed.  As  the  source  receiver  offset  chosen  for 
the  comparison  of  the  methods  was  4  kilometers  the 
asymptotic  approximation  to  the  Bessel  function  was  used. 
As  the  reflectivity  of  the  (PITL)  medium  is  frequency 
dependent  the  numerical  integration  was  carried  out  at 
543  equally  spaced  points  in  the  frequency  domain 
spanning  the  range  of  frequencies  where  the  amplitude 
spectrum  of  the  source  pulse  was  non  zero.  The  Cooley- 
Tukey  (1965)  FFT  algorithm  was  used  to  invert  to  the 
time  domain.  The  resultant  seismic  trace  is  shown  in 
Figure  2. 

The  second  trace  in  Figure  2  is  the  seismogram 
computed  via  the  ray  method.  All  possible  rays  up  to 
those  with  a  maximum  of  18  ray  segments  are  generated 
and  if  they  arrive  within  the  specified  time  window  they 
are  included  in  the  seismogram.  The  major  arrivals  and 
their  amplitudes  are  identified  along  the  time  axis. 

The  time  dependence  of  the  source  pulse  used  is 

f(t)  =  sin  (2TrfQt)  exp 


2tt  f  t 

2 

o 

y 

(4.15) 


—  co<  t  <  CO 
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Figure  4.2  Comparison  of  the  two  methods  of  computing 
synthetic  seismograms;  (a)  numerical 
integration  (b)  asymptotic  ray  theory. 
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where  fQ  is  the  predominant  frequency  of  the  pulse  and  y 
is  the  damping  factor.  In  case,  f  =  30  Hz  and  y  =  6. 
The  Fourier  transform  of  (4.  )  is 


(4.16) 


L(w)  =  - 


•  h 

ITT  y 
2tt  f 


exp  - 


2 


,  \  i 


0) 


2  TT  f 


sinh 


2 

y  a) 
47Tf 


This  pulse  was  chosen  because  of  its  simplicity  in  the 
frequency  domain;  however,  any  pulse  may  be  incorporated 
into  the  programs. 

It  is  evident  that  the  two  methods  discussed  are 
equivalent  in  this  highly  idealized  situation  in  which 
only  two  alternating  layers  are  used.  Backus  (1962)  has 
shown  that  a  medium  composed  of  many  thin  layers,  not 
necessarily  periodic,  may  be  approximated  in  the  long 
wave  length  limit  by  an  (HTI )  layer  provided  a  certain 
number  of  stability  criteria  are  satisfied. 

The  comparison  of  seismograms  of  the  above  mentioned 
type  of  thin  layered  medium  and  its  equivalent  (HTI) 
medium  is  precluded  at  the  present  time  due  to  the  large 
amounts  of  CPU  time  required.  As  an  example  for  the  case 
considered  here  the  numerical  integration  procedure 
required  250  times  the  CPU  time  as  the  ray  method,  even 
with  the  diagonalization  simplification  in  computing  the 
reflectivity. 
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Table  4.1  Specification  of  the  Media 
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APPENDIX  A 


If  the  volume  of  a  ray  tube  which  is  truncated  at 
both  ends  by  the  wavefronts  at  t  and  t  is  denoted  by  y, 
Gauss '  integral  theorem  has 


(Al) 


Hi 

V*pady  = 

j 

• 

p ( t) a • da  + 

► 

. 

Y 

Z 

(t; 

E 

t 

p  (to)a*dat 


p (t) V (t) dat~ 


p(to)V(to)dat 


I  (t) 


E(to> 


where 


da  -  incremental  surface  area  of  the  wavefront  at 
V(t)  -  normal  velocity  of  the  wavefront  at  t 

p(x)  -  density  of  the  medium  at  x 

— ^  t 

a  -  ray  velocity. 

The  surface  integration  need  not  be  carried  out  over  the 

sides  of  the  ray  tube,  as  a,  the  ray  velocity  is  always 

tangent  to  them.  The  normal  direction  to  the  surface 

of  y  is  chosen  positive  outwards  so  that  at  x=t  ,  V(t  ) 

1  o  o 

and  da^  point  in  opposite  directions  accounting  for  the 
o 

minus  sign  in  the  second  integral  in  Al . 

The  incremental  area  da^  can  be  expressed  in  terms 
of  the  ray  coordinates  (a,B,x),  (Babich  and  Buldirev 


(1972))  as 

(A2)  da^  =  J(x)dad3 


where  the  non-negative  function  J  is  linked  to  the 
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Jacobian  of  the  transformation  of  Cartesian  coordinates 
to  ray  coordinates  by  the  relation 


D  ( x  1  / x 2  ,x3) 

D (a, 3,t) 

where  a  is  the  magnitude  of  the  ray  velocity  a.  The 
function  J(t)  for  a  particular  wavefront  at  time  t  is 
equal . 


(A3) 


j  =  i 

a 


(A4 ) 


J(t)  = 


X  XX. 

a  3 


/ %x1  8x2  3x3 

where  x  =  1 -~- 

a  l  da  9a  9a 

this  it  follows  that 


and  xD  are  similar.  From 

p 


(A5) 


2  \ 

J  =  (EH-G  ) 


-► 

E  =  x  *x  , 
a  a 


H  =  V*£S  ' 


G  =  x  *xD 
a  3 


Thus  it  follows  that  since  da.  =  ^ ^  da 

o  o 


» 

{ p ( t) V (t) J ( t) -  (t  ) V(t  ) J(tQ) }dad3 
t 

(  p  (t)  V(t)  J  (t)  )  dxdad3 

j  )  J  dx  \  > 

I(t  )  t 
o  o 


( A6 ) 
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V (x ) J  (t  )  dx 


(p(t)V(t)J(t) ) V (x)dxJ (x)dadB 
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•  /• 
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and  hence 


(A7) 


V  •  pa  = 


V (x) J (t)  dx 


(p (x)V(x) J (x) ) 
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APPENDIX  B 


REFLECTION  AND  TRANSMISSION  COEFFICIENTS 


At  an  interface  between  two  elastic  media  in  welded 
contact  certain  boundary  conditions  must  be  satisfied 
when  an  SH  ray  is  incident  at  the  interface  and  produces 
a  reflected  and  transmitted  ray.  These  conditions  are 
the  continuity  of  displacement  and  shear  stress  at  the 
interface . 

Displacement: 

(Bl)  u.  +  u  =  u. 

1  r  t 

where  the  subscripted  variables  u^,  u^  and  u  refer  to 
the  displacement  components  of  the  incident,  reflected 
and  transmitted  rays  respectively. 

Shear  Stress: 


(B2 ) 


a(  ui)  +  a(ur)  =  a(ufc) 


with 


p  -  density 


B  -  as  defined  previously  in  the 
text . 

If  the  z=0  plane  is  denoted  as  the  interface  and  the 


upper  and  lower  media  are  labelled  1  and  2,  employing 
equations  (2.7),  (2.8)  and  (2.9)  from  the  text,  yields 

the  following  set  of  two  equations  in  two  unknowns. 
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( B3 ) 


-1 


u  P1B1X]_  P2b2X2J 


TR(k) 


T(k) 


-1 


L  P1B1A1 


R(k)  and  T(k)  are  the  reflection  and  transmission 
coefficients  and 


A.k2  2  V 

x  =  JJ _ i 

Aj  B.  B. 
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CO 


With  the  change  of  variable  k  =  -  q  the  system  of 


( B4 ) 


where 


(A3) 

becomes 
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The  solutions  of  (B4)  are  given  by 


( B  5 ) 


R  (q)  = 


M1W1  "  M2W2 

M1W1  +  M2W2 


(B6 ) 


T  (q)  = 


2M  w 


+  M2W2 
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Multiplying  both  the  numerator  and  denominator  of 
(B5)  by  M^w^  -  M2W2  resu-*-ts  the  following  alternate 
expression  for  the  reflection  coefficient 


(B7)  R  (q)  =  C-^iq)  -  C2(q)w2 


with 


Cx  (q) 


(M1W1 ) 2  +  (M2w2^  2 

(MjW^2  -  (M2w2)2 


C2(q) 


2M1M2w1 

(Miwi) 2  -  (M2w2)2 
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